function City_Regrs_20231117

  % t = month
  % c = area
  %
  % K  = number of coefs, including intercept if regression has one
  % N  = total number of observations across metro areas and months
  % Nt = number of observations (metro areas with data) for month t
  % Nc = number of observations (months) for area i
  %
  %                                  
  % c.not_demeaned                   Panel    SE2  = {sum c {sum t e2ct}} / (N - K)
  %                                  TS       SE2c = {sum t e2ct} / (Nc - K)                        SE2 = sum c {Nc SE2c / N}
  %                                  FM       SE2t = {sum c e2ct} / (Nt - K)                        SE2 = sum t {Nt SE2t / N}
  %								  
  % c.demean_CS_w_int or _wo_int     Panel    SE2  = {sum t [Nt/(Nt - 1)] {sum c e2ct}} / (N - K)      
  %                                  TS       SE2c = {sum t [Nt/(Nt - 1)]e2ct} / (Nc - K)           SE2 = sum c {Nc SE2c / N}
  %								  
  % c.demean_TS_w_int or _wo_int     Panel    SE2  = {sum c [Nc/(Nc - 1)] {sum t e2ct}} / (N - K)
  %                                  FM       SE2t = {sum c [Nc/(Nc - 1)]e2ct} / (Nt - K)           SE2 = sum t {Nt SE2t / N}
  %
  % c.demean_both_CS_1st or _TS_1st  Panel    SE2  = {sum c {sum t [Nc/(Nc - 1)][Nt/(Nt - 1)]e2ct}} / (N - K)
  %

  tic
   
  computer = 2;
  if computer == 0 % mac
    addpath(genpath('/users/fefama/Dropbox/Z/Papers/Real Estate'))
  elseif computer == 1 % pc
    addpath(genpath('C:\Users\fefama\Dropbox\Z\Papers\Real Estate'))
  elseif computer == 2 % pc
    addpath('d:\dropbox\projects\Real Estate', 'd:\dropbox\projects\Real Estate\matlab', 'd:\dropbox\projects\Real Estate\Data 20230103' )
  end

  c.Jan_2000 = 157;                % 200001
  c.Oct_2022 = 430;                % 202210
  
  p.last_mth = c.Oct_2022; 
  p.CPI_SA   = true;			   % true: real prices and rents  divided by CPI_SA

  c = setup_global_constants(c, p);
  p = setup_parameters(c, p);



  % setup output file...p.fid is output file
  %
  p = setup_output_file_etc(c, p);



  % load data
  %
  d = load_CPI_price_and_rent_data(c, p);


  % write raw data
  %
  if p.wrt_raw_data
    wrt_raw_data(d, c, p);
  end


  % compute and write summary stats for change in price and change in rent
  %
  cmpt_and_wrt_sum_stats(c.Jan_1987, c.last_mth, c, d, p); 


  % Run Fixed Effects, FM, Panel, and TS regressions
  %
  FE_on_FE_TS        = run_FE_on_FE_TS_regs(c, d, p);

  FM_not_DM          = run_FM_dR_on_PR_and_one_dP(c.not_demeaned,     c, d, p);
  FM_DM_CS_wo        = run_FM_dR_on_PR_and_one_dP(c.demean_CS_wo_int, c, d, p);

  Panel_not_DM       = run_Panel_dR_on_PR_and_one_dP(c.not_demeaned,     c, d, p);
  Panel_DM_CS_wo_int = run_Panel_dR_on_PR_and_one_dP(c.demean_CS_wo_int, c, d, p);

  TS_not_DM          = run_TS_dR_on_PR_and_one_dP(c.not_demeaned,        c, d, p);
  TS_DM_CS_w         = run_TS_dR_on_PR_and_one_dP(c.demean_CS_w_int,     c, d, p);

  [FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS] = add_regr_data_labels(FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS);

  control_wrt_summary_tbl_202212(FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS, c, d, p);
  
  toc
end

%-----------------------------------------------------------------------------------------
%                                                                                        *
%                                        Summary Statistics                              *
%                                                                                        *
% ****************************************************************************************
%
% Looks OK  YZ
%
% Needed
%
function cmpt_and_wrt_sum_stats(beg_mth, end_mth, c, d, p)

  
  frst_ob = [beg_mth, beg_mth + [1, 3, 6, 12], beg_mth + [1, 3, 6, 12]];
  last_ob = [end_mth, end_mth * ones(1, 4),    end_mth * ones(1, 4)];
  
 
  if p.use_nominal_prices_and_rents
    trait_data = cat(3, d.price_rent, d.chg_price, d.chg_rent);                   % c.last_mth x c.num_areas_both x (1 + 4 + 4)  
  else
    trait_data = cat(3, d.price_rent, d.chg_real_price, d.chg_real_rent);         % c.last_mth x c.num_areas_both x (1 + 4 + 4)  
  end

  TS_StdDev = squeeze(std(trait_data, 'omitnan'));
  
  Ave_TS_StdDev_all   = mean(TS_StdDev, 'omitnan');
  Ave_TS_StdDev_no_US = mean(TS_StdDev(2:end, :), 'omitnan');
  TS_StdDev_US        = TS_StdDev(1, :);
  

  fprintf(p.fid, '\n\n TS Standard Deviations of PR and Changes in Rents and Prices');

  fprintf(p.fid, '\n                    PR       dP1     dP3     dP6     dP12      dR1     dR3     dR6     dR12');
  
  fprintf(p.fid, '\nAve All Areas  %8.2f%10.2f%8.2f%8.2f%8.2f%8.2f%10.2f%8.2f%8.2f%8.2f%8.2f', Ave_TS_StdDev_all);
  
  fprintf(p.fid, '\nAve Without US %8.2f%10.2f%8.2f%8.2f%8.2f%8.2f%10.2f%8.2f%8.2f%8.2f%8.2f', Ave_TS_StdDev_no_US);
  
  fprintf(p.fid, '\nJust US        %8.2f%10.2f%8.2f%8.2f%8.2f%8.2f%10.2f%8.2f%8.2f%8.2f%8.2f', TS_StdDev_US);


  dum = NaN(1, 9);
  for trait = 1:9
    tmp = cmpt_dum_r2(squeeze(trait_data(frst_ob(trait):last_ob(trait), :, trait)), c, d, p);
    dum(trait) = tmp.r2_CS_dmnd;
  end

  fprintf(p.fid, '\n\nAdj R2 for dummy variable regressions with monthly fixed effects');

  fprintf(p.fid, '\nR2 Mthly FE    %8.2f%10.2f%8.2f%8.2f%8.2f%8.2f%10.2f%8.2f%8.2f%8.2f%8.2f', dum);

end
%
%-----------------------------------------------------------------------------------------
%                                                                                        *
%                                     End Summary Statistics                             *
%                                                                                        *
% ****************************************************************************************

% ****************************************************************************************
%                                                                                        *
%                       Functions to set up constants, load data, etc.                   *
%                                                                                        *
% ****************************************************************************************
%
% Needed
%
function c = setup_global_constants(c_in, p)

  c = c_in;
  
  c.num_2    = 2;
  c.num_3    = 3;
  c.num_4    = 4;
  c.num_5    = 5;
  c.missing  = -99.99;

  c.all      = 1;

  c.price    = 1;
  c.rent     = 2;
  c.ratio    = 3;
					 
  c.mth      = 1;
  c.qtr      = 2;
  c.hlf      = 3;
  c.ann      = 4;
  
  c.period   = [1, 3, 6, 12];
  c.num_periods = size(c.period, 2);
  
  c.nom      = 1;
  c.real     = 2;
  
  c.price_to_rent = 1;
  c.rent_to_price = 2;

  c.FM       = 1;
  c.TS       = 2;
  c.Panel    = 3;
  c.Dummy    = 4;
  c.FE_TS    = 5;
  c.FE_CS    = 6;
					
  c.slope    = 1;
  c.SE_iid   = 2;
  c.SE_LHS_1 = 3;
  c.t_iid    = 4;
  c.t_LHS_1  = 5;
  
  c.f82      = 1;
  c.f83      = 2;
  
  c.same     = 1;
  c.diff     = 2;
					
  c.header   = 1;
  c.footer   = 2;

  c.summary_1 = 1;
  c.summary_2 = 2;

  c.with_cpi = true;
  c.wout_cpi = false;
  
  c.wrt_big_t_diffs      = true;
  c.dont_wrt_big_t_diffs = false;

  c.wrt_frst_coef        = true;
  c.dont_wrt_frst_coef   = false;
  
  c.wrt_covs_etc         = true;
  c.dont_wrt_covs_etc    = false;

  c.wrt_autos            = true;
  c.dont_wrt_autos       = false;
  
  c.wrt_frst_last_mth      = true;
  c.dont_wrt_frst_last_mth = false;

  c.PR_reqd              = true;
  c.PR_not_reqd          = false;
  
  c.global_var           = true;
  c.global_stddev        = false;

  c.check_se_resid_cor      = true;
  c.dont_check_se_resid_cor = false;
  

  % Demean both x and y unless only x or only y
  %
  % Demean_CS makes sense for TS and panel regrs...set average CS value for each month to zero...can include intercept in TS
  %
  % Demean_TS makes sense for FM and panel regrs...set average TS value for each area to zero...can include intercept in FM
  %
  % Demean_both makes sense only for panel...can include overall intercept
  %        TS 1st means set average TS value for each area  to zero, then set average CS value for each month to zero (TS average no longer zero)
  %        CS 1st means set average CS value for each month to zero, then set average TS value for each area  to zero (CS average no longer zero)
  %
  % If only x or only y, include intercept
  %
  %
  % t = month
  % c = area
  %
  % K  = number of coefs, including intercept if regression has one
  % N  = total number of observations across metro areas and months
  % Nt = number of observations (metro areas with data) for month t
  % Nc = number of observations (months) for area i
  %
  %                                  
  % c.not_demeaned                   Panel    SE2  = {sum c {sum t e2ct}} / (N - K)
  %                                  TS       SE2c = {sum t e2ct} / (Nc - K)                        SE2 = sum c {Nc SE2c / N}
  %                                  FM       SE2t = {sum c e2ct} / (Nt - K)                        SE2 = sum t {Nt SE2t / N}
  %								  
  % c.demean_CS_w_int or _wo_int     Panel    SE2  = {sum t [Nt/(Nt - 1)] {sum c e2ct}} / (N - K)      
  %                                  TS       SE2c = {sum t [Nt/(Nt - 1)]e2ct} / (Nc - K)           SE2 = sum c {Nc SE2c / N}
  %								  
  % c.demean_TS_w_int or _wo_int     Panel    SE2  = {sum c [Nc/(Nc - 1)] {sum t e2ct}} / (N - K)
  %                                  FM       SE2t = {sum c [Nc/(Nc - 1)]e2ct} / (Nt - K)           SE2 = sum t {Nt SE2t / N}
  %
  % c.demean_both_CS_1st or _TS_1st  Panel    SE2  = {sum c {sum t [Nc/(Nc - 1)][Nt/(Nt - 1)]e2ct}} / (N - K)
  %
  c.not_demeaned          = 1;
  c.demean_CS_w_int       = 2;        % This makes sense only for TS and panel...include intercept
  c.demean_CS_wo_int      = 3;        % This makes sense only for TS and panel...without intercept
  c.demean_TS_w_int       = 4;        % This makes sense only for FM and panel...include intercept
  c.demean_TS_wo_int      = 5;        % This makes sense only for FM and panel...without intercept
  c.demean_both_CS_1st_w  = 6;        % This is only for panel...demean both, start with CS, then TS, include intercept
  c.demean_both_CS_1st_wo = 7;        % This is only for panel...demean both, start with CS, then TS, without intercept
  c.demean_both_TS_1st_w  = 8;        % This is only for panel...demean both, start with TS, then CS, include intercept
  c.demean_both_TS_1st_wo = 9;        % This is only for panel...demean both, start with TS, then CS, without intercept
  c.demean_both_full      = 10;       % This is only for panel...demean both, subtracts TS means and CS means...not sequential
  c.demeaned              = 11;       % This is works for dummies

  c.demean_label = ["Not demeaned", "Demean CS w int", "Demean CS wo int", "Demean TS w int", "Demean TS wo int", "Demean both, CS 1st, with int", "Demean both, CS 1st, wout int", ...
                    "Demean both, TS 1st, with int", "Demean both, TS 1st, wout int", "Demeaned both, full", "Demeaned"];

  c.CPI_Rent_Primary                 = 1;
  c.CPI_Rent_Owners_Equiv_Primary    = 2;
  c.CPI_Rent_Owners_Equiv_Residences = 3;
  
  % Assume index is observed at the end of the file month.
  %
  c.Jan_1987 = 1;     % 198701
  c.Jan_2000 = 157;   % 200001


  c.frst_mth     = c.Jan_1987;
  c.last_mth     = p.last_mth;
  c.num_mths_ttl = c.last_mth;
  c.base_date    = 198701;

  c.rescale_1987_last_mth_ave  =  0;
  c.rescale_frst_5_yr_ave      = -1;
  c.rescale_frst_10_yr_ave     = -2;


  % c.area_abrev and codes for areas (e.g., c.USAve = 1) must be linked one-for-one
  %
  %                  1       2       3        4        5         6       7         8        9       10       11       12       13       14       15      16       17       18       19
  % Include in Rent  Y       Y       Y        Y        Y         Y       Y         Y        Y       Y        Y         N       N        N         N       N        N        N        N        N        N
  c.area_abrev = ['USAve'; 'B_MA '; 'NY_NY'; 'M_FL '; 'D_MI '; 'C_IL '; 'A_GA '; 'D_TX '; 'S_WA '; 'SF_CA'; 'LA_CA'; 'SD_CA'; 'P_AZ '; 'D_CO '; 'LV_NV'; 'W_DC '; 'P_OR '; 'C_NC '; 'C_OH '; 'T_FL '; 'M_MN '];

  c.USAve = 1;  c.BMA = 2;  c.NYNY = 3;  c.MFL = 4;   c.DMI = 5;  c.CIL = 6;  c.AGA = 7;  c.DTX = 8;  c.SWA = 9;  c.SFCA = 10;  c.LACA = 11;
  c.SDCA = 12;  c.PAZ = 13; c.DCO = 14; c.LVNV = 15; c.WDC = 16; c.POR = 17; c.CNC = 18; c.COH = 19;  c.TFL = 20;  c.MMN = 21;


  % c.order_of_areas is the order the data are stored; Cities with both prices and rents MUST be first; Cities with rents but not prices are not loaded.
  %                              1        2       3      4      5      6     7      8       9      10      11
  c.order_of_areas         = [c.USAve, c.BMA, c.NYNY, c.MFL, c.DMI, c.CIL, c.AGA, c.DTX, c.SWA, c.SFCA, c.LACA];  % include only areas that are used

  c.order_for_output_price = [c.USAve, c.BMA, c.NYNY, c.MFL, c.DMI, c.CIL, c.AGA, c.DTX, c.SWA, c.SFCA, c.LACA];                                                                            % include only areas that are used
  c.order_for_output_rent  = [c.USAve, c.BMA, c.NYNY, c.MFL, c.DMI, c.CIL, c.AGA, c.DTX, c.SWA, c.SFCA, c.LACA];                                                                            % include only areas that are used
  
  c.num_areas_used    = size(c.order_of_areas, 2);
  c.num_areas_price_in = 21;
  c.num_areas_price    = 11;
  c.num_areas_rent     = 11;
  c.num_areas_both     = 11;

  c.max_TS_CS_obs      = c.num_areas_both * c.num_mths_ttl;


  % order_of_areas tells the order of areas in the data...which area comes next...from location to area...Cities with both prices and rents MUST be first
  % order_for_output tells the order of areas in the output...which area comes next...from location to area
  % area_loc_in_data tells where a area is in the data...which column...from area to location
  % area_loc_in_output tells where a area is in the output...which row...from area to location

  c.area_loc_in_data = zeros(1, c.num_areas_used);
  for area = 1:c.num_areas_used
    c.area_loc_in_data(c.order_of_areas(area)) = area;  % where is a area in the data...which column for a area
  end
  
  % c.full_area_name must line up with c.area_abrev and codes for areas (e.g., c.USAve = 1) must be linked one-for-one
  %
  c.full_area_name = strings(c.num_areas_price, 1);
  
  c.full_area_name(1,  1) = "US Average";
  c.full_area_name(2,  1) = "Boston";
  c.full_area_name(3,  1) = "New York";
  c.full_area_name(4,  1) = "Miami";
  c.full_area_name(5,  1) = "Detroit";
  c.full_area_name(6,  1) = "Chicago";
  c.full_area_name(7,  1) = "Atlanta";
  c.full_area_name(8,  1) = "Dallas";
  c.full_area_name(9,  1) = "Seattle";
  c.full_area_name(10, 1) = "San Francisco";
  c.full_area_name(11, 1) = "Los Angeles";

  c.full_area_name_char(1,  :) = 'US Average   ';
  c.full_area_name_char(2,  :) = 'Boston       ';
  c.full_area_name_char(3,  :) = 'New York     ';
  c.full_area_name_char(4,  :) = 'Miami        ';
  c.full_area_name_char(5,  :) = 'Detroit      ';
  c.full_area_name_char(6,  :) = 'Chicago      ';
  c.full_area_name_char(7,  :) = 'Atlanta      ';
  c.full_area_name_char(8,  :) = 'Dallas       ';
  c.full_area_name_char(9,  :) = 'Seattle      ';
  c.full_area_name_char(10, :) = 'San Francisco';
  c.full_area_name_char(11, :) = 'Los Angeles  ';

end

% Needed
%
function p = setup_parameters(c, p)

  p.cc_changes                         = true;                      % this affects all price changes in program; cc or percent changes, eg ln[p(t+1)/P(t)] or [p(t+1)-P(t)]/P(t)
  p.use_nominal_prices_and_rents       = false;
  
  p.SA_prices                          = true;                      % seasonally adjusted prices or not
  p.SA_rents                           = false;                     % seasonally adjusted rents or not; Relevant only for All
  p.PR_or_RP                           = c.price_to_rent;           % c.price_to_rent or c.rent_to_price
  
  p.PR_lag                             = 0;                       % Lag for PR in regressions
  
  p.subtract_grnd_mn_both_full         = true;                     % Subtract overall mean when demeaning both_full
  
  p.P_and_R_reqd_for_all_mthly_means   = false;                    % If true, include only areas with both price and rent when computing all monthly means
  p.P_and_R_reqd_for_all_area_means    = false;                    % If true, include only areas with both price and rent when computing all monthly means
  
  p.sum_stats_frst_ob_max              = true;                      % If true, frst_ob is max of frst price, frst_rent, and beg; if false, frst_ob differs for price and rent
  
  p.summary_includes_FM_Std            = true;
  p.summary_includes_Panel_TS_CS_1st   = false;
  
  % at least one of these must be true
  %
  p.incld_same_beg_end                 = false;                     % if true, include set of regressions with the same start and same finish month; all use the latest beg and earliest end in the set
  p.incld_diff_beg_end                 = true;                      % if true, include set of regressions that each uses the earliest start and latest finish as possible

  p.frst_last_form_mth_same_or_diff    = c.diff;
  
  p.ave_chg_req_price_and_rent         = true;

  p.out_format                         = c.f83;                     % c.f82 or c.f83...%8.3f
  p.slope_format                       = c.f82;                     % c.f82 or c.f83...%8.3f
  p.se_format                          = c.f82;                     % c.f82 or c.f83...%8.3f
  p.t_stat_format                      = c.f82;                     % c.f82 or c.f83...%8.3f

  if ~p.incld_same_beg_end && ~p.incld_diff_beg_end
    disp('either p.incld_same_beg_end or p.incld_diff_beg_end (or both) must be true')
    pause
  end
  
  p.LHS_RHS_vars                       = [c.rent,  c.price; ...  % dependent variable
                                          c.ratio, c.ratio; ...  % 1st independent variable...used in univariate regression
                                          c.price, c.price];     % 2nd independent variable...used in all but first regression in set
  
  
  p.global_var_or_stddev               = c.global_var;

  if p.incld_same_beg_end
    p.frst_same_diff = c.same;
  else
    p.frst_same_diff = c.diff;
  end

  if p.incld_diff_beg_end
    p.last_same_diff = c.diff;
  else
    p.last_same_diff = c.same;
  end
  
  p.skip_auto_adj_if_sum_neg = true;       % if true, use iid SE if adj SE for auto_cors or cross_cors would be smaller
  
  
  % parameters for FM regressions
  %
  % rescale date matters only for P/R and only in FM regressions
  % rescaling has no effect on pcnt or cc change in rents and prices, and rescaling is absorbed in mean of P/R (into intercept) in TS regression
  % rescaling does affect relative values of P/R in FM regressions; 
  % it controls for differences in P and R indices across areas caused by differences in scaling of P and R and by differences in quality, tax rates, etc
  %
  p.rescale_date                      = c.rescale_1987_last_mth_ave;  %NaN;  %c.Jan_2000; % c.rescale_1987_last_mth_ave; % c.Jan_2000;         % NaN means don't rescale      % p.rescale_date must be a specific month or c.rescale_1987_last_mth_ave
  p.rescale_date_txt                  = "Ave";  %"None"; % "Jan 2000"; % "Ave";

  p.wrt_limited_output                = true;
  p.wrt_TS_area_output                = true;

  p.wrt_FM_dR_on_PR_and_one_dP_coefs  = false;

  p.wrt_same_frst_last_form_mths      = false;

  p.wrt_SE_adj                        = true;
  
  p.num_top_t_diffs                   = 10;
  
  p.FM_DM_both_inclds_int             = false;
  

  %
  p.wrt_raw_data                          = false;                    % this is for diagnostics
  
  p.num_coef_stats_simple                 = 5;

  % parameters for FM regressions of dR on P/R and dP
  %
  % There are three formats
  %   1. Similar to the new format TS results .........................................called "new" here too
  %   2. Single dP or no dP (0, 1, 3, 6, or 12) with period determined by dP...........called "one_diff"
  %   3. Single dP or no dP (0, 1, 3, 6, or 12) with period determined by longest dP...called "one_same"
  %
  %
  p.ln_PR_in_regrs	                      = true;   % false;
  p.wrt_FM_dR_on_PR_and_one_dP_autos      = false;  %true;

  % General parameters
  %
  p.LHS_per                               = [1, 3, 6, 12];
  p.RHS_per                               = [1, 3, 6, 12];
  
  p.summary_LHS_pers                      = [3, 12];
  p.summary_LHS_per_loc                   = [2, 4];

  p.summary_RHS_pers                      = [3, 12];
  p.summary_RHS_per_loc                   = [2, 4];


  % General FM parameters
  %
  p.FM_max_coefs_allowed_in_regr          = 5;                % this is the maximum number of explanatory variables in any FM cross-section regression; some regressions have only seven areas
  

  % General Panel parameters
  %
  p.panel_time_clstrd_just_contemp        = false;

  p.summary_Panel_Both_wo_int             = true;     % include intercept in Panel_Both regressions in summary at bottom of output

  % General TS parameters
  %
  p.num_autos_resids                      = 36;
  p.num_autos_FM_coefs                    = 36;        % must be at least 1
  p.num_autos_sum_stats                   = 36;
  p.num_autos_in_FM_coef_se               = 0;         %[12, 24, 36];
  p.num_autos_in_TS_coef_se               = 0;         %[12, 24, 36];


  % TS "one" format is dR on PR and one dP
  %
  p.wrt_TS_dR_on_PR_and_one_dP_autos      = false;
  p.wrt_TS_dR_on_PR_and_one_dP_summary    = true;                    % This is abbreviated output at the bottom
  
  p.num_autos_FE_resids                   = 12;

  p.TS_DM_both_inclds_int                 = false;

end
  
% Needed
%
function p = setup_output_file_etc(c, p)

  body = "FF_City_Regrs_20231117_";
  
  if p.use_nominal_prices_and_rents
    nominal_or_real = "nominal";
  else
    nominal_or_real = "real";
  end


  label_last = "Oct_2022";
  
  rescale_txt = "PR_rescaled_by_ave_1987_2022110";
  
  if p.ln_PR_in_regrs  
    ln_PR_txt = "ln_PR";
  else
    ln_PR_txt = "not_ln_PR";
  end
  
  PR_lag_txt = sprintf('PR_lag_%d', p.PR_lag);

  if p.summary_Panel_Both_wo_int
    sum_Pnl_Both_wo_txt = "Pnl_Bth_wo";
  else 
    sum_Pnl_Both_wo_txt = "Pnl_Bth_w";
  end

  
  if p.subtract_grnd_mn_both_full
    grand_mean = "zero_grnd_mn_both_full";
  else
    grand_mean = "ignore_grnd_mn_both_full";
  end
  
  filename = sprintf('%s%s_%s_%s_%s_%s_%s_%s_test.txt', body, nominal_or_real, label_last, ln_PR_txt, PR_lag_txt, rescale_txt, sum_Pnl_Both_wo_txt, grand_mean);
  
  p.fid = fopen(filename, 'wt');
  
  fprintf(p.fid, '\nOutput produced by City_Regrs_202301030.m');

  if p.SA_prices
    fprintf(p.fid, '\n\nPrices are seasonally adjusted.');
  else
    fprintf(p.fid, '\n\nPrices are not seasonally adjusted.');
  end
  
  if p.SA_rents
    fprintf(p.fid, '\nRents for All is seasonally adjusted; Seasonally adjusted data are not available for individual areas.');
  else
    fprintf(p.fid, '\nRents for All and individual areas are not seasonally adjusted.');
  end
  
  if p.cc_changes
    fprintf(p.fid, '\n\nChanges in prices and rents are continuously compounded.');
  else	
    fprintf(p.fid, '\n\nChanges in prices and rents are simple, not continuously compounded.');
  end
  
  if p.use_nominal_prices_and_rents
    fprintf(p.fid, '\n\nThe regressions use nominal changes in prices and rents.');
  else
    fprintf(p.fid, '\n\nThe regressions use real (inflation adjusted) changes in prices and rents.');
    if p.CPI_SA
      fprintf(p.fid, '\n\nWe use SEASONALLY ADJUSTED CPI-U to convert from nominal to real.');
    else
      fprintf(p.fid, '\n\nWe use NON SEASONALLY ADJUSTED CPI-U to convert from nominal to real.');
    end
  end
 	
  fprintf(p.fid, '\n\nInflation rates and changes in prices and rents are multiplied by 100.');
  
 
  if p.PR_or_RP == c.price_to_rent
    fprintf(p.fid, '\n\nThis run uses price/rent, not rent/price.');
  else
    fprintf(p.fid, '\n\nThis run uses rent/price, not price/rent.');
  end
  
  if p.ln_PR_in_regrs
    fprintf(p.fid, '\n\nThe ratios of price to rent (or rent to price) in the regressions are in logs.');
  else
    fprintf(p.fid, '\n\nThe ratios of price to rent (or rent to price) in the regressions are NOT in logs.');
  end    
  
  if p.PR_lag == 0
    fprintf(p.fid, '\n\nPR is not lagged in the regressions.');
  elseif p.PR_lag == 1
    fprintf(p.fid, '\n\nPR is lagged one month in the regressions. NOTE...We do not lag PR in the summary statistics.');
  else
    fprintf(p.fid, '\n\nPR is lagged %d months in the regressions. NOTE...We do not lag PR in the summary statistics.', p.PR_lag);
  end
  
  fprintf(p.fid, '\n\nFM regressions can have no more than %d explanatory variables.', p.FM_max_coefs_allowed_in_regr);
  
  fprintf(p.fid, '\n\nRescale date matters only for PR (or RP) and only in FM regressions. Rescaling has no effect on percent');
  fprintf(p.fid, '\nor cc changes in rents and prices, and rescaling of PR is absorbed into the intercept in TS regressions.');
  fprintf(p.fid, '\nRescaling does affect relative values of PR in FM regressions. We set each area''s average value of PR');
  fprintf(p.fid, '\nfor the full 1987-2022 sample to 1. The goal is to control for cross-area differences in PR caused by');
  fprintf(p.fid, '\ndifferences in scaling of P and R and by permanent differences in quality, tax rates, etc.');
  fprintf(p.fid, '\n\nThe program can rescale to a specific month, so all values of  PR or R/P are 1 in that month. We cannot run FM');
  fprintf(p.fid, '\nregressions that include PR in the rescale month. The missing month complicates calculation of autocorrelations.');
  fprintf(p.fid, '\nAs a work-around when we rescale to a specific month, we use the averages of the other months as the coefficients');
  fprintf(p.fid, '\nfor the missing month.');
  
  fprintf(p.fid, '\n\nPrices, rents, price/rent, and rent/price are rescaled so averages for 1987 to 202210 are 1.0');

end
  
% Needed
%
function d = load_CPI_price_and_rent_data(c, p)


  % Read CPI for January 1987 to Nov 2022
  %
  CPI = dlmread('CPI_Monthly_194701_202211.txt', '', [481, 0, 481 + c.last_mth - 1, 2]);   % date,  NSA,  SA

  if p.CPI_SA     % true: real prices and rents  divided by CPI_SA
    d.CPI = CPI(:, 3);
  else
    d.CPI = CPI(:, 2);
  end

  d.idate   = CPI(:, 1);


  % setup price and rent data
  %
    % get input file names
  %
  [price_filename, rent_filename] = get_input_file_names(c, p);

  d.price     = NaN(c.num_mths_ttl, c.num_areas_price);
  d.rent      = NaN(c.num_mths_ttl, c.num_areas_both);
  
  d.area_name  = reshape(blanks(c.num_areas_used * 13), c.num_areas_used, 13);
  d.area_abrev = reshape(blanks(c.num_areas_used * 5),  c.num_areas_used, 5);

  d.frst_ob   = zeros(c.num_areas_used, 3);
  d.last_ob   = zeros(c.num_areas_used, 3);

  % read price data  
  %
  f_in = fopen(price_filename, 'r');
  
  area_name_in  = textscan(f_in, [repmat('%s',1,c.num_areas_price_in), '%*[^\n]'], 1, 'HeaderLines', 34); 
  area_abrev_in = textscan(f_in, [repmat('%s',1,c.num_areas_price_in), '%*[^\n]'], 1); 
  frst_ob_in   = textscan(f_in, [repmat('%d',1,c.num_areas_price_in), '%*[^\n]'], 1);
  last_ob_in   = textscan(f_in, [repmat('%d',1,c.num_areas_price_in), '%*[^\n]'], 1);

  area_name  = reshape(blanks(c.num_areas_price_in * 13), c.num_areas_price_in, 13);
  area_abrev = reshape(blanks(c.num_areas_price_in * 5),  c.num_areas_price_in, 5);
  frst_ob   = cell2mat(frst_ob_in);
  last_ob   = cell2mat(last_ob_in);
  
  
  for row = 1:c.num_areas_price_in
    area_name (row, :) = strjust(sprintf('%13s', cell2mat(area_name_in {1, row})), 'left');
    area_abrev(row, :) = strjust(sprintf('%5s',  cell2mat(area_abrev_in{1, row})), 'left');
    frst_ob(row)       = loc_mth(frst_ob(row), c);
    last_ob(row)       = loc_mth(last_ob(row), c);
  end

  
  price_in = textscan(f_in, ['%d', repmat('%f',1,c.num_areas_price_in), '%*[^\n]'], c.num_mths_ttl, 'CollectOutput',1);
  price_idate = price_in{1};
  price       = price_in{2};

  loc_has_data = zeros(c.num_areas_used, 1);  

  for area = 1:c.num_areas_price_in
    loc = find_area_abrev_match(area_abrev(area, :), c);

    if loc == 0
      disp('Cant find area abreviation');
      disp(area);
      disp(area_abrev(area, :));
      pause;

    elseif loc <= c.num_areas_both
      loc_has_data(c.area_loc_in_data(loc)) = true;
      
      d.price     (:, c.area_loc_in_data(loc))       = price    (:, area);    % c.area_loc_in_data is where the area is in c.order of areas
      d.frst_ob   (c.area_loc_in_data(loc), c.price) = frst_ob  (area);       % frst_ob(:, 1) and last_ob(:, 1) is for price
      d.last_ob   (c.area_loc_in_data(loc), c.price) = min(last_ob(area), c.last_mth) ;
      d.area_name (c.area_loc_in_data(loc), :)       = area_name (area, :);
      d.area_abrev(c.area_loc_in_data(loc), :)       = area_abrev(area, :);
    end
  end
  fclose(f_in);

  % read rent data  
  %
  f_in = fopen(rent_filename, 'r');
  
  area_name_in  = textscan(f_in, [repmat('%s',1,c.num_areas_rent + 1), '%*[^\n]'], 1, 'HeaderLines', 1); 
  area_abrev_in = textscan(f_in, [repmat('%s',1,c.num_areas_rent + 1), '%*[^\n]'], 1); 
  frst_ob_in   = textscan(f_in, [repmat('%d',1,c.num_areas_rent + 1), '%*[^\n]'], 1);
  last_ob_in   = textscan(f_in, [repmat('%d',1,c.num_areas_rent + 1), '%*[^\n]'], 1);

  area_abrev = reshape(blanks(c.num_areas_rent * 5), c.num_areas_rent, 5);
  frst_ob   = cell2mat(frst_ob_in);
  last_ob   = cell2mat(last_ob_in);
  
  for row = 1:c.num_areas_rent + 1
    area_name (row, :) = strjust(sprintf('%13s', cell2mat(area_name_in {1, row})), 'left');
    area_abrev(row, :) = strjust(sprintf('%5s',  cell2mat(area_abrev_in{1, row})), 'left');
    frst_ob  (row)    = loc_mth(frst_ob(row), c);
    last_ob  (row)    = min([loc_mth(last_ob(row), c), c.last_mth]);
  end
  
  rent_in    = textscan(f_in, ['%d', repmat('%f',1,c.num_areas_rent + 1), '%*[^\n]'], c.num_mths_ttl, 'CollectOutput',1);
  rent_idate = rent_in{1};
  rent       = rent_in{2};
  

  Check_Alignment(d.idate, price_idate, rent_idate);
 
 
  % Load rents in columns defined by c.order_of_areas
  %
  for area = 1:c.num_areas_rent + 1
  
    if (p.SA_rents && strcmp(area_abrev(area, :),'SAAve')) || (~p.SA_rents && strcmp(area_abrev(area, :),'NSAAv'))
      loc = find_area_abrev_match('USAve', c);  % loc is the location of the area in c.area_name, which is also the value of the area's name, e.g. c.USAve = 1
      
    elseif ~strcmp(area_abrev(area, :),'SAAve') && ~strcmp(area_abrev(area, :),'NSAAv')    % want only one of SA Ave and NSA Ave; that one is selected in if statement above; skip the other
      loc = find_area_abrev_match(area_abrev(area, :), c);
    end

    if loc == 0
      disp('Cant find area abreviation');
      disp(area);
      disp(area_abrev(area, :));
      pause;

    elseif loc <= c.num_areas_both
      d.rent(:, c.area_loc_in_data(loc))         = rent(:, area);
      d.frst_ob(c.area_loc_in_data(loc), c.rent) = frst_ob(area);    % frst_ob(:, 1) and last_ob(:, 1) is for price
      d.last_ob(c.area_loc_in_data(loc), c.rent) = last_ob(area);

      if ~loc_has_data(c.area_loc_in_data(loc))  % loc_has_data = t if price data available for area
        disp('Area in rent file is in c.order_of_areas but was not found in price file.')
        disp(area_abrev(area, :))
        pause;
        
        % Here because area has no price data
        %
        loc_has_data(c.area_loc_in_data(loc))    = true;
        d.area_name (c.area_loc_in_data(loc), :) = area_name (area, :);
        d.area_abrev(c.area_loc_in_data(loc), :) = area_abrev(area, :);

      end
    end
  end
  fclose(f_in);


  % build header with area abbreviations
  %
  d.header_areas = reshape(blanks(c.num_areas_used * 8), c.num_areas_used, 8);
  for row = 1:c.num_areas_used
    d.header_areas(row, :) = sprintf('   %5s', strjust(d.area_abrev(row, :), 'right'));
  end      


  % Compute ratio of rent to price and price to rent, and the average values of the ratios (used in the const_x regression)
  %
  d.rent_price           = NaN(c.last_mth, c.num_areas_both);
  d.price_rent           = NaN(c.last_mth, c.num_areas_both);
  d.mean_rent_price_area =          NaN(1, c.num_areas_both);
  d.mean_price_rent_area =          NaN(1, c.num_areas_both);

  d.mean_rent_price_mth  = NaN(c.last_mth, 1);
  d.mean_price_rent_mth  = NaN(c.last_mth, 1);
  
  for area = 1:c.num_areas_both
    if 0 < min(d.frst_ob(area, c.price:c.rent))
      d.frst_ob(area, c.ratio) = max(d.frst_ob(area, c.price:c.rent));
      d.last_ob(area, c.ratio) = min(d.last_ob(area, c.price:c.rent));
      
      frst_mth = d.frst_ob(area, c.ratio);
      last_mth = d.last_ob(area, c.ratio);
      if frst_mth <= last_mth
      
        if p.ln_PR_in_regrs
          d.rent_price(frst_mth:last_mth, area) = log(d.rent (frst_mth:last_mth, area) ./ d.price(frst_mth:last_mth, area));
          d.price_rent(frst_mth:last_mth, area) = log(d.price(frst_mth:last_mth, area) ./ d.rent (frst_mth:last_mth, area));
        else
          d.rent_price(frst_mth:last_mth, area) = d.rent (frst_mth:last_mth, area) ./ d.price(frst_mth:last_mth, area);
          d.price_rent(frst_mth:last_mth, area) = d.price(frst_mth:last_mth, area) ./ d.rent (frst_mth:last_mth, area);
        end
                
        d.mean_price_rent_area(1, area) = mean(d.price_rent(frst_mth:last_mth, area), 'omitnan');
        d.mean_rent_price_area(1, area) = mean(d.rent_price(frst_mth:last_mth, area), 'omitnan');
      end
    end
  end    


  % Rescale PR so full period average for each area is 1
  %
  % When plotting, we rescale prices and rents so they are 1 on Jan 2000; 
  % Here the base is usually the average value from 1987-2020. This choice affects only the FM regressions. 
  % 
  for area = 1:c.num_areas_both
    
    if d.frst_ob(area, c.price) ~= 0
      area_mean_price = mean(d.price(d.frst_ob(area, c.price):d.last_ob(area, c.price), area));
      d.price(d.frst_ob(area, c.price):d.last_ob(area, c.price), area) = d.price(d.frst_ob(area, c.price):d.last_ob(area, c.price), area) ./ area_mean_price;
    end
    
    if d.frst_ob(area, c.rent) ~= 0
      area_mean_rent = mean(d.rent(d.frst_ob(area, c.rent):d.last_ob(area, c.rent), area));
      d.rent(d.frst_ob(area, c.rent):d.last_ob(area, c.rent), area) = d.rent(d.frst_ob(area, c.rent):d.last_ob(area, c.rent), area) ./ area_mean_rent;
    end

    if d.frst_ob(area, c.ratio) ~= 0
      frst_mth = d.frst_ob(area, c.ratio);
      last_mth = d.last_ob(area, c.ratio);
      if frst_mth <= last_mth
        if p.ln_PR_in_regrs
          d.price_rent(frst_mth:last_mth, area) = d.price_rent(frst_mth:last_mth, area) -  d.mean_price_rent_area(1, area);
          d.rent_price(frst_mth:last_mth, area) = d.rent_price(frst_mth:last_mth, area) -  d.mean_rent_price_area(1, area);
        else
          d.price_rent(frst_mth:last_mth, area) = d.price_rent(frst_mth:last_mth, area) ./ d.mean_price_rent_area(1, area);
          d.rent_price(frst_mth:last_mth, area) = d.rent_price(frst_mth:last_mth, area) ./ d.mean_rent_price_area(1, area);
        end
      end
    end
  end

  d.CPI = d.CPI ./ mean(d.CPI(c.Jan_1987:end));
  

  % compute real prices and rents
  %
  d.real_price = d.price ./ d.CPI;
  d.real_rent  = d.rent  ./ d.CPI;


  % Compute growth rates
  %
  d.chg_price      = NaN(c.last_mth, 4, c.num_areas_both);
  d.chg_rent       = NaN(c.last_mth, 4, c.num_areas_both);
  d.chg_real_price = NaN(c.last_mth, 4, c.num_areas_both);
  d.chg_real_rent  = NaN(c.last_mth, 4, c.num_areas_both);

  for area = 1:c.num_areas_both
    if 0 < d.frst_ob(area, c.rent)
           d.chg_rent(:, :, area) = control_chg_x(d.rent     (:, area), d.frst_ob(area, c.rent), d.last_ob(area, c.rent), c, p);
      d.chg_real_rent(:, :, area) = control_chg_x(d.real_rent(:, area), d.frst_ob(area, c.rent), d.last_ob(area, c.rent), c, p);
    end

    if 0 < d.frst_ob(area, c.price)   
           d.chg_price(:, :, area) = control_chg_x(d.price     (:, area), d.frst_ob(area, c.price), d.last_ob(area, c.price), c, p);
      d.chg_real_price(:, :, area) = control_chg_x(d.real_price(:, area), d.frst_ob(area, c.price), d.last_ob(area, c.price), c, p);
    end
  end


  % Change order to mth, area, period
  %
  d.chg_rent       = permute(d.chg_rent,       [1, 3, 2]);    
  d.chg_real_rent  = permute(d.chg_real_rent,  [1, 3, 2]);
  d.chg_price      = permute(d.chg_price,      [1, 3, 2]);
  d.chg_real_price = permute(d.chg_real_price, [1, 3, 2]);


  % Average and SDs of changes in price and rent are for periods used in regressions.
  % Changes are dated in the month they are observed.
  %
  % Mth:  Price and rent changes are available from d.frst_ob(area, c.price + 1) to d.last_ob(area, c.price) and d.frst_ob(area, c.rent + 1) to d.last_ob(area, c.price)
  %
  %       Regression for month t is R(t+1)/R(t) on PR(t) and P(t-1)/P(t), so t can go from Max[d.frst_ob(area, c.price) + 1, d.frst_ob(area, c.rent)] to Min[d.last_ob(area, c.price), d.last_ob(area, c.rent) - 1]
  %
  %       PR is available for every t with P(t)/P(t-1) and R(t+1)/R(t)
  %
  % Qtr:  Regression for month t is R(t+3)/R(t) on PR(t) and P(t-3)/P(t), so t can go from Max[d.frst_ob(area, c.price) + 3, d.frst_ob(area, c.rent)] to Min[d.last_ob(area, c.price), d.last_ob(area, c.rent) - 3]
  %
  [d.CS_ave_chg_price,      d.CS_SD_chg_price,      d.CS_n_chg_price] = compute_mth_means_and_SDs(d.chg_price,      d.chg_price, d.chg_rent, c, d, p);                    % output is c.last_mth x c.num_areas_both x 4
  [d.CS_ave_chg_real_price, d.CS_SD_chg_real_price, d.CS_n_chg_price] = compute_mth_means_and_SDs(d.chg_real_price, d.chg_price, d.chg_rent, c, d, p);                    % output is c.last_mth x c.num_areas_both x 4
						   
  [d.CS_ave_chg_rent,       d.CS_SD_chg_rent,       d.CS_n_chg_rent]  = compute_mth_means_and_SDs(d.chg_rent,       d.chg_rent,  d.chg_rent, c, d, p);                    % output is c.last_mth x c.num_areas_both x 4
  [d.CS_ave_chg_real_rent,  d.CS_SD_chg_real_rent,  d.CS_n_chg_rent]  = compute_mth_means_and_SDs(d.chg_real_rent,  d.chg_rent,  d.chg_rent, c, d, p);                    % output is c.last_mth x c.num_areas_both x 4
						   
						   
  [d.TS_ave_chg_price,      d.TS_SD_chg_price,      d.TS_n_chg_price] = compute_area_means_and_SDs_full(d.chg_price,      d.chg_price, d.chg_rent, c.rent,  c, d, p); % output is c.num_areas_both x 4
  [d.TS_ave_chg_real_price, d.TS_SD_chg_real_price, d.TS_n_chg_price] = compute_area_means_and_SDs_full(d.chg_real_price, d.chg_price, d.chg_rent, c.rent,  c, d, p); % output is c.num_areas_both x 4
						   
  [d.TS_ave_chg_rent,       d.TS_SD_chg_rent,       d.TS_n_chg_rent]  = compute_area_means_and_SDs_full(d.chg_rent,       d.chg_rent,  d.chg_rent, c.price, c, d, p); % output is c.num_areas_both x 4
  [d.TS_ave_chg_real_rent,  d.TS_SD_chg_real_rent,  d.TS_n_chg_rent]  = compute_area_means_and_SDs_full(d.chg_real_rent,  d.chg_rent,  d.chg_rent, c.price, c, d, p); % output is c.num_areas_both x 4

  [d.TS_ave_rent_price, d.TS_SD_rent_price, d.TS_n_rent_price] = compute_area_means_and_SDs_full(d.rent_price, d.chg_price, d.chg_rent, c.ratio, c, d, p);            % output is c.num_areas_both x 4  d.chg_price_mth and d.chg_rent_mth are ignored for c.ratio
  [d.TS_ave_price_rent, d.TS_SD_price_rent, d.TS_n_price_rent] = compute_area_means_and_SDs_full(d.price_rent, d.chg_price, d.chg_rent, c.ratio, c, d, p);            % output is c.num_areas_both x 4

  d.CS_ave_rent_price = mean(d.rent_price, 2, 'omitnan');
  d.CS_ave_price_rent = mean(d.price_rent, 2, 'omitnan');
  d.CS_SD_rent_price  =  std(d.rent_price, 0, 2, 'omitnan');
  d.CS_SD_price_rent  =  std(d.price_rent, 0, 2, 'omitnan');
  d.CS_n_price_rent   = sum(~isnan(d.price_rent), 2);
  
end

% Needed
%  
function chg_all = control_chg_x(x, frst_ob, last_ob, c, p)
  
  chg_all = NaN(c.last_mth, 4);
  
  if 0 < frst_ob
    gap = [1, 3, 6, 12];
   
    for gap_indx = 1:size(gap, 2)
      chg_all(:, gap_indx) = chg_x(x, frst_ob, last_ob, gap(gap_indx), c, p);
    end
  end
  
end

% Needed
%
function output = chg_x(input, frst_ob, last_ob, gap, c, p)

  output = NaN(c.last_mth, 1);
  
  if p.cc_changes
    output(frst_ob + gap:last_ob, 1) = 100 * log(input(frst_ob + gap:last_ob, 1) ./ input(frst_ob:last_ob - gap, 1));
  else
    output(frst_ob + gap:last_ob, 1) = 100 *   ((input(frst_ob + gap:last_ob, 1) ./ input(frst_ob:last_ob - gap, 1)) - 1);
  end
end   

  
function [mthly_means, mthly_SDs, mthly_cnts] = compute_mth_means_and_SDs(x, x_price, x_rent, c, d, p)    % output is c.num_mths_ttl x 1

  y = x;
  
  if p.P_and_R_reqd_for_all_mthly_means
    y((isnan(x_price) + isnan(x_rent)) > 0) = NaN;
%%    y((isnan(x_price(:, 1:c.num_areas_both, :)) + isnan(x_rent)) > 0) = NaN;
  end

  mthly_means = mean(y, 2, 'omitnan');
  mthly_SDs   =  std(y, 0, 2, 'omitnan');
  mthly_cnts  = sum(~isnan(y), 2);
end  


function [area_means, area_SDs, area_cnts] = compute_area_means_and_SDs_full(x, x_price, x_rent, price_or_rent, c, d, p)     % output is 1 x c.num_areas_both x 1 or 4

  % Average and SDs of changes in price and rent are for periods used in regressions.
  % Changes are dated in the month they are observed.
  %
  % Mth:  Price and rent changes are available from d.frst_ob(area, c.price) + 1 to d.last_ob(area, c.price) and d.frst_ob(area, c.rent + 1) to d.last_ob(area, c.price)
  %
  %       Regression for month t is R(t+1)/R(t) on PR(t) and P(t-1)/P(t), so t can go from Max[d.frst_ob(area, c.price) + 1, d.frst_ob(area, c.rent)] to Min[d.last_ob(area, c.price), d.last_ob(area, c.rent) - 1]
  %
  %       PR is available for every t with P(t)/P(t-1) and R(t+1)/R(t)
  %
  % Qtr:  Regression for month t is R(t+3)/R(t) on PR(t) and P(t-3)/P(t), so t can go from Max[d.frst_ob(area, c.price) + 3, d.frst_ob(area, c.rent)] to Min[d.last_ob(area, c.price), d.last_ob(area, c.rent) - 3]
  %
  y = x;
  
  % price or rent
  %
  if price_or_rent == c.ratio
    num_periods = 1;
    
  else
    num_periods = c.num_periods;
    
    if p.P_and_R_reqd_for_all_area_means
      y((isnan(x_price(:, 1:c.num_areas_both, :)) + isnan(x_rent)) > 0) = NaN;
    end
  end

  area_means =   NaN(1, size(y, 2), num_periods);
  area_SDs   =   NaN(1, size(y, 2), num_periods);
  area_cnts  = zeros(1, size(y, 2), num_periods);

  for mth_qtr_etc = 1:num_periods

    for area = 1:c.num_areas_both
      frst_mth = d.frst_ob(area, price_or_rent) + c.period(mth_qtr_etc); 
      last_mth = d.last_ob(area, price_or_rent); 

      area_means(1, area, mth_qtr_etc) =       mean(y(frst_mth:last_mth, area, mth_qtr_etc), 'omitnan');
        area_SDs(1, area, mth_qtr_etc) =        std(y(frst_mth:last_mth, area, mth_qtr_etc), 'omitnan');
       area_cnts(1, area, mth_qtr_etc) = sum(~isnan(y(frst_mth:last_mth, area, mth_qtr_etc)));
    end
  end
  
end  
  

% Needed
%
function [price_filename, rent_filename] = get_input_file_names(c, p)

  if p.SA_prices
    price_filename = 'CS_Prices_SA_198701_202210.txt';
  else
    price_filename = 'CS NSA Prices 198701_202210.txt';
    disp('Non-seasonal CS Prices were not downloaded.')
    pause
  end
  
  rent_filename = 'CPI_Rent_of_Primary_Residence_NSA_198701_202211.txt';

end
  
% ****************************************************************************************
%                                                                                        *
%                     End Functions to Set Up Constants, Load Data, etc.                 *
%                                                                                        *
% ****************************************************************************************
%
% ****************************************************************************************
%                                                                                        *
%                       Utility Functions used in Old and/or New                         *
%                                                                                        *
% ****************************************************************************************
%
% Needed
%
function loc = loc_mth(date, c)   % date is yyyymm, like 201012   loc is the location relative to the c.base_date, like 198701
  
  cal_yr   = floor(date ./ 100);
  cal_mth  = mod(date, 100);
  base_yr  = floor(c.base_date ./ 100);
  base_mth = mod(c.base_date, 100);

  loc      = 12 * (cal_yr - base_yr) + (cal_mth - base_mth) + 1;
end
  
% Needed
%
function loc = find_area_abrev_match(area_abrev, c)
  
  for loc = 1:size(c.area_abrev, 1)
    if strcmp(c.area_abrev(loc, :), area_abrev)
      return
    end
  end
  loc = 0;
end

% Needed
%
function wrt_raw_data(d, c, p)

  fprintf(p.fid, '\n\n\nCS Prices      (wrt_raw_data)');
  fprintf(p.fid, '\n\n      '); fprintf(p.fid, '  %10s', d.header_areas');

  for j = 1:size(d.price, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%10.4f', d.price(j, :));
  end
  
  fprintf(p.fid, '\n\n\nCPI Rents and CPI            (wrt_raw_data)');
  fprintf(p.fid, '\n\n      ');   fprintf(p.fid, '  %10s', d.header_areas(1:c.num_areas_both, :)'); fprintf(p.fid, '       CPI');
  
  for j = 1:size(d.rent, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%10.4f', d.rent(j, 1:c.num_areas_both), d.CPI(j));
  end

  fprintf(p.fid, '\n\n\nPrice/Rent            (wrt_raw_data)');
  fprintf(p.fid, '\n\n      ');   fprintf(p.fid, '  %10s', d.header_areas(1:c.num_areas_both, :)'); fprintf(p.fid, '       CPI');
  
  for j = 1:size(d.rent, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%10.4f', d.price_rent(j, 1:c.num_areas_both), d.CPI(j));
  end

  fprintf(p.fid, '\n\n\nRent/Price            (wrt_raw_data)');
  fprintf(p.fid, '\n\n      ');   fprintf(p.fid, '  %10s', d.header_areas(1:c.num_areas_both, :)'); fprintf(p.fid, '       CPI');
  
  for j = 1:size(d.rent, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%10.4f', d.rent_price(j, 1:c.num_areas_both), d.CPI(j));
  end

end

% Needed
%
function Check_Alignment(idate, price_idate, rent_idate)

  % Make sure files are aligned
  %
  if size(rent_idate, 1) ~= size(price_idate, 1)
    disp('sizes of rent_idate and price_idate do not match')
    size(rent_idate), size(price_idate)
    disp('rent 1 & end...price 1 & end')
    [rent_idate(1), rent_idate(end), price_idate(1), price_idate(end)]
    pause;
  end
  
  if size(rent_idate, 1) ~= size(idate, 1)
    disp('sizes of rent_idate and cpi idate do not match')
    size(rent_idate), size(idate)
    disp('rent 1 & end...idate 1 & end')
    [rent_idate(1), rent_idate(end), idate(1), idate(end)]
    pause;
  end
  
  for j = 1:size(idate, 1)
    if rent_idate(j, 1) ~= price_idate(j, 1)
      disp('rent_idate and price_idate misaligned')
      rent_idate(j, 1), price_idate(j, 1)
      pause;
    elseif idate(j, 1) ~= rent_idate(j, 1)
      disp('rent_idate and d.idate misaligned')
      rent_idate(j, 1), idate(j, 1)
      pause;
    end
  end

end  


%-----------------------------------------------------------------------------------------
%                                                                                        *
%                       Utility Functions for Running Regressions                        *
%                                                                                        *
% ****************************************************************************************
%
% Needed
%
function chg = load_mthly_qtrly_semi_or_annual_chg(price_or_rent, per, area, c, d, p)

  % use area = 0 for all
  %
  indx_per = find(c.period == per);
  
  if p.use_nominal_prices_and_rents && price_or_rent == c.price
    chg = d.chg_price(:, :, indx_per);
  elseif p.use_nominal_prices_and_rents && price_or_rent == c.rent
    chg = d.chg_rent(:, :, indx_per);
  elseif ~p.use_nominal_prices_and_rents && price_or_rent == c.price
    chg = d.chg_real_price(:, :, indx_per);
  elseif ~p.use_nominal_prices_and_rents && price_or_rent == c.rent
    chg = d.chg_real_rent(:, :, indx_per);
  end

  if area ~= 0
    chg = chg(:, area);
  end
end
  
function x = load_full_PR_vector_regr(area, c, d, p)

  if area == 0
    if p.PR_or_RP == c.rent_to_price
      x = NaN(size(d.rent_price));
      x(p.PR_lag + 1:end, :) = d.rent_price(1:end - p.PR_lag, :);
    else
      x = NaN(size(d.price_rent));
      x(p.PR_lag + 1:end, :) = d.price_rent(1:end - p.PR_lag, :);
    end
  else
    if p.PR_or_RP == c.rent_to_price
      x = NaN(size(d.rent_price(:, area)));
      x(p.PR_lag + 1:end, 1) = d.rent_price(1:end - p.PR_lag, area);
    else    
      x = NaN(size(d.price_rent(:, area)));
      x(p.PR_lag + 1:end, 1) = d.price_rent(1:end - p.PR_lag, area);
    end
  end
end


function num_regrs = count_TS_regrs_for_ttl_RHS_per(ttl_RHS_per, p)

  num_regrs = 1;
  
  for RHS_per_indx = 1:size(p.RHS_per, 2)
    RHS_per = p.RHS_per(RHS_per_indx);
          
    if (RHS_per <= ttl_RHS_per) && (mod(ttl_RHS_per, RHS_per) == 0)
      num_regrs = num_regrs + 1;
    end
  end
end



% OK  Dates are correct, nt and nc are correct
%
function [frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, nt, nc] = determine_and_wrt_frst_and_last_form_mth(LHS_per, LHS_var, RHS_var, wrt_frst_last_mth_passed, c, d, p)

  % Computes first and last formation month for each area.
  %
  %   if p.frst_last_form_mth_same_or_diff == c.diff...frst_form_mth_area(c.num_areas, 1:5), last_form_mth_area(c.num_areas, 1)  frst_form_mth(1:5, 1) = min(frst_form_mth_area)  last_form_mth = max(last_form_mth_area)
  %
  %   if p.frst_last_form_mth_same_or_diff == c.same...frst_form_mth_area(c.num_areas, 1),   last_form_mth_area(c.num_areas, 1)  frst_form_mth         = min(frst_form_mth_area)  last_form_mth = max(last_form_mth_area)
  %      RHS_per is max(p.RHS_per)
  %   

  RHS_per = [0, p.RHS_per];
  
  num_RHS_pers  = size(p.RHS_per, 2);
  max_RHS_per   =  max(p.RHS_per);
  
  frst_form_mth = NaN(1 + num_RHS_pers, 1); 
  last_form_mth = NaN;                      

  frst_form_mth_area = NaN(c.num_areas_both, 1 + num_RHS_pers);    % 1 + num_RHS_pers...1 is for simple P/R
  last_form_mth_area = NaN(c.num_areas_both, 1);                   
  
  nt = zeros(c.num_mths_ttl,   num_RHS_pers + 1);                  % how many areas for month
  nc = zeros(c.num_areas_both, num_RHS_pers + 1);                  % how many months for area
  
  % Compute frst and last formation month for each area
  %

  if p.frst_last_form_mth_same_or_diff == c.diff

    [frst_form_mth(:), last_form_mth] = cmpt_frst_and_last_formation_mths(0, LHS_var, RHS_var, LHS_per, RHS_per, c, d, p, c.PR_reqd);  % frst_form_mth 1 x num_RHS_pers, last_form_mth scalar determined by LHS_per

    for area = 1:c.num_areas_both
      [frst_form_mth_area(area, :), last_form_mth_area(area, 1)] = cmpt_frst_and_last_formation_mths(area, LHS_var, RHS_var, LHS_per, RHS_per, c, d, p, c.PR_reqd);  % frst_form_mth 1 x (1 + num_RHS_pers), last_form_mth scalar determined by LHS_per
       frst_form_mth_area(area, :) = max([squeeze(frst_form_mth_area(area, :)); frst_form_mth(:)']);  % note: frst_form_mth_same_area can vary across RHS periods
       last_form_mth_area(area, 1) = min([last_form_mth_area(area, 1); last_form_mth]);   % last_form_mth_diff is actually scalar; last_form_mth_diff_area(area, 1) should always be <= last_form_mth_diff
    end
    
  else

    [frst_form_mth(1), last_form_mth] = cmpt_frst_and_last_formation_mths(0, LHS_var, RHS_var, LHS_per, max_RHS_per, c, d, p, c.PR_reqd);  % frst_form_mth scalar, last_form_mth scalar determined by LHS_per
    frst_form_mth(:) = frst_form_mth(1);
  
    for area = 1:c.num_areas_both
      frst_form_mth_area(area, :) = max(max(frst_form_mth_area(area, :)), frst_form_mth(1));   % note: frst_form_mth_same_area can vary across RHS periods
      last_form_mth_area(area, 1) = min([last_form_mth_area(area, 1); last_form_mth]);      % last_form_mth_diff_area(area, 1) should always be <= last_form_mth_diff
    end
  end
  
  for area = 1:c.num_areas_both
    for RHS_per_indx = 1:1 + num_RHS_pers
      nc(area, RHS_per_indx) = last_form_mth_area(area) - frst_form_mth_area(area, RHS_per_indx) + 1;
      
      nt(frst_form_mth_area(area, RHS_per_indx):last_form_mth_area(area), RHS_per_indx) = nt(frst_form_mth_area(area, RHS_per_indx):last_form_mth_area(area), RHS_per_indx) + 1;
    end
  end

  % write frst and last formation months for each area
  %
  if wrt_frst_last_mth_passed

    fprintf(p.fid, '\n\nFirst and last formation months for each area');
    if p.frst_last_form_mth_same_or_diff == c.diff
      fprintf(p.fid, '\n                        Diff'); 
    else 
      fprintf(p.fid, '\n                        Same'); 
    end
    fprintf(p.fid, '\n          0       1       3       6      12     Last');

    for area = 1:c.num_areas_both
      fprintf(p.fid, '\n%5s', c.area_abrev(c.order_of_areas(area), :)); fprintf(p.fid, '%8d', d.idate(frst_form_mth_area(area, :)), d.idate(last_form_mth_area(area, 1)));
    end
  end
end

% frst_form_mth is a function of RHS_per, and RHS_var   1 x num_RHS_pers   usually 1 + size(p.RHS_per)
% last_form_mth is a function of LHS_per and LHS_var    scalaar
%
function [frst_form_mth, last_form_mth] = cmpt_frst_and_last_formation_mths(area, LHS_var, RHS_var, LHS_per, RHS_per, c, d, p, PR_reqd_passed)
%
% If LHS_var = c.ratio, LHS_per and RHS_per are irrelevant
%
% Note: frst_form_mth is 1 x num_RHS_pers and last_form_mth is 1 x num_LHS_pers; must test TS_min_mths <= last_form_mth - frst_form_mth + 1 in application
%
% Logic: frst_form_mth must be late  enough to have good RHS and LHS changes and, if with_PR_passed, to compute ratio
%        last_form_mth must be early enough to have good RHS and LHS changes and, if with_PR_passed, to compute ratio
%
% If PR_reqd_passed
%  frst_form_mth = Max(d.frst_ob(area, RHS_var) + RHS_per, d.frst_ob(area, LHS_var))               % ensures RHS and LHS data don't start after first month of either RHS and LHS change, combination of two conditions ensures good ratio
%  last_form_mth = Min(d.last_ob(area, RHS_var),           d.last_ob(area, LHS_var) - LHS_per)     % ensures RHS and LHS data don't stop before last  month of either RHS and LHS change, combination of two conditions ensures good ratio
% else
%  frst_form_mth = Max(d.frst_ob(area, RHS_var) + RHS_per, d.frst_ob(area, LHS_var))               % ensures RHS and LHS data don't start after first month of either RHS and LHS change
%  last_form_mth = Max(d.last_ob(area, RHS_var),           d.last_ob(area, LHS_var) - LHS_per)		 % ensures RHS and LHS data don't stop before last  month of either RHS and LHS change
%
% RHS_per and LHS_per are 1 x num_RHS_pers and 1 x num_LHS_pers, so frst_RHS_mth is 1 x num_RHS_pers and last_LHS_mth is 1 x num_LHS_pers.
%
% Note: frst_form_mth is 1 x num_RHS_pers and last_form_mth is 1 x num_LHS_pers; must test TS_min_mths <= last_form_mth - frst_form_mth + 1 in application
%
  if LHS_var == 0
    disp("It makes no sense to have LHS_var = 0 in cmpt_frst_and_last_formation_mths")
    pause
  end
  
  if RHS_var == 0 && ~PR_reqd_passed
    disp("It makes no sense to have both RHS_var == 0 and PR_reqd_passed == false")
    pause
  end

  if area == 0

    if RHS_var ~= 0
      frst_form_mth = min(d.frst_ob(1:c.num_areas_both, RHS_var)) + RHS_per;
    else
      frst_form_mth = 0;
    end
    
    last_form_mth = max(d.last_ob(1:c.num_areas_both, LHS_var)) - LHS_per;
  
    % get first and last month with both P and R for any area
    %
    if PR_reqd_passed
      frst_form_mth = max(frst_form_mth, min(max(d.frst_ob(1:c.num_areas_both, :)')) + p.PR_lag);
      last_form_mth = min(last_form_mth, max(min(d.last_ob(1:c.num_areas_both, :)')) + p.PR_lag);
    end
    
    
  else
    last_form_mth = d.last_ob(area, LHS_var) - LHS_per;

    if RHS_var == 0
      frst_form_mth = min(d.frst_ob(area, :));
    else
      frst_form_mth = d.frst_ob(area, RHS_var) + RHS_per;
      last_form_mth(last_form_mth > d.last_ob(area, RHS_var)) = d.last_ob(area, RHS_var);         % RHS data can't stop before last formation month
    end
    
    frst_form_mth(frst_form_mth < d.frst_ob(area, LHS_var)) = d.frst_ob(area, LHS_var);   % LHS data can't start more than gap months after first formation mth..if this dominates, frst_form_mth can be constant for c.diff

    if PR_reqd_passed                                                                       % must have price and rent in every formation month
      frst_form_mth(frst_form_mth < max(d.frst_ob(area, :)) + p.PR_lag) = max(d.frst_ob(area, :) + p.PR_lag);
      last_form_mth(last_form_mth > min(d.last_ob(area, :)) + p.PR_lag) = min(d.last_ob(area, :) + p.PR_lag);
    end
    
  end
end

function [frst_form_mth, last_form_mth] = cmpt_frst_and_last_formation_mths_just_dummies(area_in, LHS_var, LHS_per, c, d)

  RHS_var = mod(LHS_var, 2) + 1;
  
  if area_in == 0
  
    % find first form month for any area and last form month for any area
    %
    frst_form_mth = c.last_mth;
    last_form_mth = 0;
    for area = 1:c.num_areas_both
      frst_form_mth = min(frst_form_mth, max(d.frst_ob(area, :))+ p.PR_lag);                                        % this is the first month area has both price and rent
      last_form_mth = max(last_form_mth, min(d.last_ob(area, :) + p.PR_lag, d.last_ob(area, LHS_var) - LHS_per));   % this is the earlier of last month area has both price and rent and last formation month with LHS change available
    end

  else  
    frst_form_mth = max(d.frst_ob(area_in, :) + p.PR_lag);                                                      % this is the first month area has both price and rent
    last_form_mth = min(d.last_ob(area_in, :) + p.PR_lag, d.last_ob(area_in, LHS_var) - LHS_per);               % this is the earlier of last month area has both price and rent and last formation month with LHS change available
  end
end

% ****************************************************************************************
%                                                                                        *
%                       Shared Functions for Regression Output                           *
%                                                                                        *
% ****************************************************************************************
function wrt_header_for_regr_results(FM_Panel_or_TS, LHS_var, RHS_var, style_passed, c, d, p)  

  if FM_Panel_or_TS == c.FM
    FM_Panel_or_TS_label = 'FM';
  elseif FM_Panel_or_TS == c.Panel
    FM_Panel_or_TS_label = 'Panel';
  else
    FM_Panel_or_TS_label = 'TS';
  end

  PR_label       = load_PR_label(c, p);
  LHS_full_label = load_P_or_R_label(LHS_var, c);
  RHS_full_label = load_P_or_R_label(RHS_var, c);
  

  fprintf(p.fid, '\n\n\n%s regrs forecasting change in %s with %s, a past change in %s, and a future change in %s    (wrt_header_for_regr_results, %s, dR_on_PR_and_one_dP)', FM_Panel_or_TS_label, LHS_full_label, PR_label, RHS_full_label, RHS_full_label, FM_Panel_or_TS_label);

                
  if FM_Panel_or_TS == c.FM
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nThe variables are not demeaned and each monthly regression includes an intercept.');
    elseif style_passed == c.demean_TS_w_int
      fprintf(p.fid, '\nEach area''s variables are demeaned and each monthly regression includes an intercept.');
    elseif style_passed == c.demean_TS_wo_int
      fprintf(p.fid, '\nEach area''s variables are demeaned and the monthly regressions do not include an intercept.');
    end
  
    describe_adj_t_stats_FM(p);
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nAve var e is the average of each month''s average squared residual, with each residual adjusted for the degrees of freedom used estimating its month''s coefficients.');
      fprintf(p.fid, '\nGlobal var e is the average squared residual, with each residual adjusted for the degrees of freedom used estimating its month''s coefficients.');
    else
      fprintf(p.fid, '\nAve var e is the average of each month''s average squared residual, with each residual adjusted first for demeaning then for the degrees of freedom used estimating its month''s coefficients.');
      fprintf(p.fid, '\n        Adj e2 = e2 * [n_area / (n_area - 1)] / (n_mth - n_coefs)');
      fprintf(p.fid, '\nGlobal var e is the global average of the squared adjusted residuals.');
    end
    fprintf(p.fid, '\nThe R2 for month t = 1 - month t''s average adjusted squared residual / var of month t''s raw Y.');
    fprintf(p.fid, '\nThe Global R2 is 1 - average of all months'' adjusted squared residuals / var of all raw Y.');
    
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nave_SD_x is the TS average of the monthly CS standard deviations of the raw xs.');
      fprintf(p.fid, '\nSD_SD_x is the TS standard deviation of the monthly CS standard deviations of the raw xs.');
    else
      fprintf(p.fid, '\nave_SD_x is the TS average of the monthly CS standard deviations of the demeaned xs.');
      fprintf(p.fid, '\nSD_SD_x is the TS standard deviation of the monthly CS standard deviations of the demeaned xs.');
    end
    
    
  elseif FM_Panel_or_TS == c.Panel
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nThe variables are not demeaned and the panel regression includes an intercept.');
    elseif style_passed == c.demean_CS_w_int
      fprintf(p.fid, '\nEach monthly CS of variables is demeaned and the panel regression includes an intercept.');
    elseif style_passed == c.demean_CS_wo_int
      fprintf(p.fid, '\nEach monthly CS of variables is demeaned and the panel regression does not include an intercept.');
    elseif style_passed == c.demean_TS_w_int
      fprintf(p.fid, '\nEach area TS of variables is demeaned and the panel regression includes an intercept.');
    elseif style_passed == c.demean_TS_wo_int
      fprintf(p.fid, '\nEach area TS of variables is demeaned and the panel regression does not include an intercept.');
    
    elseif style_passed == c.demean_both_CS_1st_w
      fprintf(p.fid, '\nAll the area TS and monthly CS of variables are demeaned, the cross-sections are demeaned first, and the panel regression includes an intercept.');
    elseif style_passed == c.demean_both_CS_1st_wo
      fprintf(p.fid, '\nAll the area TS and monthly CS of variables are demeaned, the cross-sections are demeaned first, and the panel regression does not include an intercept.');
    elseif style_passed == c.demean_both_TS_1st_w
      fprintf(p.fid, '\nAll the area TS and monthly CS of variables are demeaned, the time-series are demeaned first, and the panel regression includes an intercept.');
    elseif style_passed == c.demean_both_TS_1st_wo
      fprintf(p.fid, '\nAll the area TS and monthly CS of variables are demeaned, the time-series are demeaned first, and the panel regression does not include an intercept.');
    end
    
  
    describe_adj_t_stats_Panel(p);
    
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nGlobal Res Var is the over-all sum of squared residuals divided by N-K (the total number of observations minus the number of coefficients in the regression).');

    elseif x_in_y(style_passed, [c.demean_CS_w_int, c.demean_CS_wo_int])
      fprintf(p.fid, '\nGlobal Res Var is the global average of the adjusted squared residuals, with each squared residual adjusted for monthly CS');
      fprintf(p.fid, '\ndemeaning then for the degrees of freedom used estimating the panel coefficients.');
      fprintf(p.fid, '\n        Adj Res2 = Res2 * [n_mth / (n_mth - 1)] * [N / (N - K)');

    elseif x_in_y(style_passed, [c.demean_TS_w_int, c.demean_TS_wo_int])
      fprintf(p.fid, '\nGlobal Res Var is the global average of the adjusted squared residuals, with each squared residual adjusted for area TS');
      fprintf(p.fid, '\ndemeaning then for the degrees of freedom used estimating the panel coefficients.');
      fprintf(p.fid, '\n        Adj Res2 = Res2 * [n_area / (n_area - 1)] * [N / (N - K)');

    elseif x_in_y(style_passed, [c.demean_both_CS_1st_w, c.demean_both_CS_1st_wo, c.demean_both_TS_1st_w, c.demean_both_TS_1st_wo])
      fprintf(p.fid, '\nGlobal var e is the global average of the adjusted squared residuals, with each squared residual adjusted for monthly CS');
      fprintf(p.fid, '\ndemeaning, then for area TS demeaning, then for the degrees of freedom used estimating the panel coefficients.');
      fprintf(p.fid, '\n        Adj Res2 = Res2 * [n_mth / (n_mth - 1)] * [n_area / (n_area - 1)] * [N / (N - K)');
    end
    fprintf(p.fid, '\nGlobal R2 is 1 - Global Res Var / Var estimated with all raw Y.');
    
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nave_SD_x is the TS average of the monthly CS standard deviations of the raw xs. (The Panel and FM regression definitions match.)');
      fprintf(p.fid, '\nSD_SD_x is the TS standard deviation of the monthly CS standard deviations of the raw xs.');
    else
      fprintf(p.fid, '\nave_SD_x is the TS average of the monthly CS standard deviations computed after each x is demeaned by its area''s average. (The Panel and FM regression definitions match.)');
      fprintf(p.fid, '\nSD_SD_x is the TS standard deviation of the monthly CS standard deviations computed after each x is demeaned by its area''s average.');
    end

  elseif FM_Panel_or_TS == c.TS
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nThe variables are not demeaned and each area regression includes an intercept.');
    elseif style_passed == c.demean_CS_w_int
      fprintf(p.fid, '\nEach month''s variables are demeaned and each area''s regression includes an intercept.');
    elseif style_passed == c.demean_CS_wo_int
      fprintf(p.fid, '\nEach month''s variables are demeaned and the area regressions do not include an intercept.');
    end
  
    describe_adj_t_stats_TS(p);
    
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nAve var e is the cross-area average of the areas'' TS average squared residuals, with each residual adjusted for the degrees of freedom');
      fprintf(p.fid, '\nused estimating its area''s coefficients.');
      fprintf(p.fid, '\nGlobal var e is the global average of the adjusted squared residuals, with each squared residual adjusted for the degrees of freedom used');
      fprintf(p.fid, '\nestimating its area''s coefficients.');
    else
      fprintf(p.fid, '\nAve var e is the cross-area average of the areas'' TS average adjusted squared residuals, with each residual adjusted first for monthly');
      fprintf(p.fid, '\ndemeaning then for the degrees of freedom used estimating its area''s coefficients.');
      fprintf(p.fid, '\n        Adj e2 = e2 * [n_mth / (n_mth - 1)] *[n_area / (n_area - n_coefs)]');
      fprintf(p.fid, '\nGlobal var e is the global average of the squared adjusted residuals.');
    end
    fprintf(p.fid, '\nThe R2 for area i = 1 - area i''s average adjusted squared residual / var of area i''s raw Y.');
    fprintf(p.fid, '\nThe Global R2 is 1 - average of all areas'' adjusted squared residuals / var of all raw Y.');
    
    if style_passed == c.not_demeaned
      fprintf(p.fid, '\nave_SD_x is the CS average of the TS standard deviations of each area''s raw xs.');
      fprintf(p.fid, '\nSD_SD_x is the CS standard deviation of the TS standard deviations of each area''s raw xs.');
    else
      fprintf(p.fid, '\nave_SD_x is the CS average of the TS standard deviations of each area''s xs after each monthly CS is demeaned.');
      fprintf(p.fid, '\nSD_SD_x is the CS standard deviation of the TS standard deviations of each area''s xs after each monthly CS is demeaned.');
    end
  end 
 
  fprintf(p.fid, '\n');
    
end

function PR_label = load_PR_label(c, p)

  if p.ln_PR_in_regrs && p.PR_or_RP == c.price_to_rent
    PR_label = "ln(price/rent)";
        
  elseif p.ln_PR_in_regrs && p.PR_or_RP == c.rent_to_price
    PR_label = "ln(rent/price)";
    
  elseif ~p.ln_PR_in_regrs && p.PR_or_RP == c.price_to_rent
    PR_label = "price/rent";

  elseif ~p.ln_PR_in_regrs && p.PR_or_RP == c.rent_to_price
    PR_label = "rent/price";
  end
end
  
function PR_label_shrt = load_PR_label_shrt(c, p)
  if p.PR_or_RP == c.price_to_rent
    PR_label_shrt = "PR ";
  else
    PR_label_shrt = "RP ";
  end
end

function label = load_P_or_R_label(var_value, c)

  if var_value == c.price
    label = "price";
  elseif var_value == c.rent
    label = "rent";
  end
end


function label = load_dP_or_dR_label(var_value, c)

  if var_value == c.price
    label = " dP";
  elseif var_value == c.rent
    label = " dR";
  end
end


function wrt_one_block_Dum(label, Dum, LHS_RHS_var_indx, c, d, p)

  fprintf(p.fid, '\n\n%s', label);
  for LHS_per_indx = 1:size(p.LHS_per, 2)
    fprintf(p.fid, '\n%2d    ', p.LHS_per(LHS_per_indx)); 
    
    if p.global_var_or_stddev == c.global_stddev
      fprintf(p.fid, '%8.3f                                                                   %8.3f                                                                   %8.3f', ...
          Dum.glbl_r2_ttl_raw(LHS_per_indx, LHS_RHS_var_indx), sqrt(Dum.glbl_var_e(LHS_per_indx, LHS_RHS_var_indx)), sqrt(Dum.glbl_var_y(LHS_per_indx, LHS_RHS_var_indx)));
    else
      fprintf(p.fid, '%8.3f                                                                   %8.3f                                                                   %8.3f', ...
          Dum.glbl_r2_ttl_raw(LHS_per_indx, LHS_RHS_var_indx), Dum.glbl_var_e(LHS_per_indx, LHS_RHS_var_indx), Dum.glbl_var_y(LHS_per_indx, LHS_RHS_var_indx));
    end
  end
end

function wrt_one_block(label, Panel, LHS_RHS_var_indx, c, d, p)

  if size(Panel.glbl_r2_ttl_raw, 1) == 5
    blnk = '                                   ';
  else
    blnk = '   ';
  end

  fprintf(p.fid, '\n\n%s', label);
  for LHS_per_indx = 1:size(p.LHS_per, 2)
    fprintf(p.fid, '\n%2d    ', p.LHS_per(LHS_per_indx)); 
    fprintf(p.fid, '%8.3f', Panel.glbl_r2_ttl_raw(:, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '%s', blnk);

    if p.global_var_or_stddev == c.global_stddev
      fprintf(p.fid, '%8.3f', sqrt(Panel.glbl_var_e    (:, LHS_per_indx, LHS_RHS_var_indx))); fprintf(p.fid, '%s', blnk);
      fprintf(p.fid, '%8.3f', sqrt(Panel.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx))); 
    else
      fprintf(p.fid, '%8.3f', Panel.glbl_var_e    (:, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '%s', blnk);
      fprintf(p.fid, '%8.3f', Panel.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx)); 
    end
  end
end


function wrt_global_r2_SD_header(LHS_RHS_var_indx, c, p)

  LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
  RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);

  PR_label_shrt = load_PR_label_shrt(c, p);
  LHS_dZ_label  = load_dP_or_dR_label(LHS_var, c);
  RHS_dZ_label  = load_dP_or_dR_label(RHS_var, c);

  fprintf(p.fid, '\n\nRegressing %s on %s and %s   (wrt_one_block)', LHS_dZ_label, PR_label_shrt, RHS_dZ_label);
    
  if p.global_var_or_stddev == c.global_stddev
    fprintf(p.fid, '\n                                Global Bias-Adjusted R2                                           Global Std Dev of Bias-Adjusted Residuals                                           Global Std Dev of Raw Y');
  else
    fprintf(p.fid, '\n                                Global Bias-Adjusted R2                                           Global Variance of Bias-Adjusted Residuals                                          Global Variance of Raw Y');
  end
    
  % fprintf(p.fid, '\nRHS                    1         3         6         12                     1         3         6         12                     1         3         6         12 ');
  % fprintf(p.fid, '\nLHS          PR      PR&dP     PR&dP     PR&dP     PR&dP          PR      PR&dP     PR&dP     PR&dP     PR&dP          PR      PR&dP     PR&dP     PR&dP     PR&dP');
  %                    1    1234567890123456789012345678901234567890123456789012312345678901234567890123456789012345678901234567890
  
  
  fprintf(p.fid, '\nRHS  '); 
  for blk = 1:3 
    fprintf(p.fid, '            %3d     %3d     %3d      %3d                                   ', p.RHS_per);                     %d         %d         %d         %d', p.RHS_per, p.RHS_per);
  end
  
  fprintf(p.fid, '\nLHS   '); 
  for blk = 1:3
     fprintf(p.fid, '     %s ', strtrim(PR_label_shrt)); 
     for regr = 1:size(p.RHS_per, 2)
       fprintf(p.fid, '   %s&%s', strtrim(PR_label_shrt), strtrim(RHS_dZ_label));
     end
     fprintf(p.fid, '    dP%2d    dP%2d    dP%2d    dP%2d   ', p.RHS_per);
   end

end                                                                              

%-----------------------------------------------------------------------------------------
%                                                                                        *
%                                 Statistical Functions                                  *
%                                                                                        *
% ****************************************************************************************
%
function [coef_stats, autos] = summarize_FM_coefs(coef, LHS_per, c, d, p)   % (3 + 2 + num_addl_ts) x num_coefs = (ave coef, se_smpl, se_LHS_1, t_smpl, t_LHS_1, addl_ts... num_autos x num_coefs

  [num_coef_stats, num_addl_ts, num_autos] = get_num_FM_coef_stats(p);   % num_coef_stats = 5 + num_addl_ts = coef, SE simpl, SE LHS-1, t-simple, t-LHS-1, addl ts (e.g., 12, 24, 36)

  num_obs   = size(coef, 1);
  num_coefs = size(coef, 2);
  t         = NaN(2 + num_addl_ts,  num_coefs);        % simple (t_0)...t_P (dP Per or RHS Per)...t_12 (1st element in p.num_autos_in_FM_coef_se)...t_24 etc...t_m 
  
  coef_ave  = mean(coef, 'omitnan');                     % 1 x num_coefs
  coef_sd   = std (coef, 'omitnan');                     % 1 x num_coefs
  se_smpl   = coef_sd ./ sqrt(num_obs);       % 1 x num_coefs
  se_LHS_1  = se_smpl;
  t(1, :)   = coef_ave ./ se_smpl;            % 1 x num_coefs
  t(2, :)   = t(1, :);

  autos = NaN(num_autos, num_coefs);

  for coef_indx = 1:num_coefs
    autos_shrt = autocorr(coef(:, coef_indx), 'NumLags', num_autos);
    autos(1:num_autos, coef_indx) = autos_shrt(2:num_autos + 1);

    % sum autos up to LHS_per - 1 
    %
    if 1 < LHS_per
      num_lags = LHS_per - 1;
      auto_adj = sqrt(cmpt_auto_adj(num_lags, num_obs, autos_shrt(2:num_lags + 1)));
      
      if 1 <= auto_adj
        se_LHS_1(coef_indx) = se_smpl(coef_indx) * auto_adj;
            t(2, coef_indx) =    t(1, coef_indx) / auto_adj;
      end
    end
    
    if num_addl_ts ~=0
      for auto_indx = 1:size(p.num_autos_in_FM_coef_se, 2)
        num_lags = p.num_autos_in_FM_coef_se(auto_indx);
        auto_adj = sqrt(cmpt_auto_adj(num_lags, num_obs, autos_shrt(2:num_lags + 1)));
        if 1 <= auto_adj
          t(2 + auto_indx, coef_indx) = t(1, coef_indx) / sqrt(cmpt_auto_adj(num_lags, num_obs, autos_shrt(2:num_lags + 1)));
        end
      end
    end
  end
  
  % (3 + num_ts + num_autos) x num_coefs = (ave coef, se_smpl, se_LHS_1, t smpl, t adj, t autos) x num_coefs
  % 
  % num_ts    = 2 + size(p.num_autos_in_FM_coef_se, 2)...simple, LHS-1, 
  % num_autos = max([LHS_per - 1, p.num_autos_in_FM_coef_se]);
  %
  coef_stats = [coef_ave; se_smpl; se_LHS_1; t];  

end  
    
function [coef_stats, ave_resid_cor, se_resid_cor_good] = summarize_TS_area_coefs(coef, resids, c, d, p)

  num_obs    = size(coef,   1);
  num_coefs  = size(coef,   2);
  num_areas = size(resids, 2);

  % compute sum of correlations of contemporaneous residuals
  %
  sum_resid_cors = sum(corrcoef(resids, 'Rows', 'pairwise'), 'all') - num_areas;  % Sum of all off-diagonal residual correlations...'Rows', 'pairwise' means skip observations if either variable does not have data
  
  ave_resid_cor  = sum_resid_cors / (num_areas * (num_areas - 1));
  
  
  coef_ave     = mean(coef);                                          % 1 x num_coefs
  coef_sd      = std (coef);                                          % 1 x num_coefs
  se_smpl      = coef_sd ./ sqrt(num_obs);                            % 1 x num_coefs
  t_smpl       = coef_ave ./ se_smpl;                                 % 1 x num_coefs

  if (p.skip_auto_adj_if_sum_neg && sum_resid_cors < 0) || sum_resid_cors / num_areas <= 0
    se_resid_cor      = se_smpl;                                          % 1 x num_coefs
    t_resid_cor       = t_smpl;                                           % 1 x num_coefs
    se_resid_cor_good = false;
  else
    se_resid_cor = se_smpl .* sqrt(1 + sum_resid_cors / num_areas);  % 1 x num_coefs
    t_resid_cor  = coef_ave ./ se_resid_cor;                          % 1 x num_coefs
    se_resid_cor_good = true;
  end
  
  coef_stats = [coef_ave; se_smpl; se_resid_cor; t_smpl; t_resid_cor];

end  

function coef_stats = summarize_TS_area_coefs_simple(coef, c, d, p)

  num_obs   = size(coef, 1);
  num_coefs = size(coef, 2);

  coef_ave   = mean(coef);                       % 1 x num_coefs
  coef_sd    = std (coef);                       % 1 x num_coefs
  se_smpl    = coef_sd ./ sqrt(num_obs);         % 1 x num_coefs
  t_smpl     = coef_ave ./ se_smpl;              % 1 x num_coefs

  coef_stats = [coef_ave; se_smpl; t_smpl];
end  

function auto_adj = cmpt_auto_adj(num_lags, num_obs, autos)

  sum_auto = 1;
  for lag = 1:num_lags
    sum_auto = sum_auto + 2 * (1 - lag ./ num_obs) .* autos(lag);
  end
  
  auto_adj = sum_auto;
end


%-----------------------------------------------------------------------------------------
%                                                                                        *
%                       Utility Functions for Writing Output                             *
%                                                                                        *
% ****************************************************************************************
%
function describe_adj_t_stats_FM(p)

  fprintf(p.fid, '\n\nThe t-statistics for the coefficients are adjusted for serial correlation.');
  fprintf(p.fid, '\nt simple is the simple t-statistic, with no adjustment for autocorrelation.');
  fprintf(p.fid, '\nt LHS-1 adjusts for autocorrelations up to the LHS period - 1. If LHS_var is the 3-month change in price, t LHS-1 uses auocorrelations up to lag 2.');
  
  num_t_lag = numel(p.num_autos_in_FM_coef_se);
  if num_t_lag == 0
  elseif num_t_lag == 1
    fprintf(p.fid, '\nt nn adjusts for autocorrelations up to lag %d', p.num_autos_in_FM_coef_se);
  elseif num_t_lag == 2
    fprintf(p.fid, '\nt nn adjusts for autocorrelations up to lags %d and %d.', p.num_autos_in_FM_coef_se);
  else
    fprintf(p.fid, '\nt nn adjusts for autocorrelations up to lags ');
    for nn = p.num_autos_in_FM_coef_se(1:end - 1); fprintf(p.fid, ' %d,', nn); end; fprintf(p.fid, ' and %d', p.num_autos_in_FM_coef_se(end));
  end
end


function describe_adj_t_stats_TS(p)

  fprintf(p.fid, '\n\nThe t-statistics for the coefficients are adjusted for serial correlation, heteroscedasticity, or both.');
  fprintf(p.fid, '\nt simple is the simple t-statistic.');
  fprintf(p.fid, '\nt LHS-1 adjusts for residual autocorrelation up to the LHS period - 1. If LHS_var is the 3-month change in price, t LHS-1 uses auocorrelations up to lag 2.');
  fprintf(p.fid, '\nt hetero allows for heteroscedasticity in the residual covariance matrix.');
  fprintf(p.fid, '\nthac L-1 allows for heteroscedasticity and residual autocorrelation up to the LHS period - 1.');
end


function describe_adj_t_stats_Panel(p)

  fprintf(p.fid, '\n\nThe t-statistics for the coefficients are adjusted for heteroscedasticity, general cross-correlation, general serial-correlation, or both.');
  fprintf(p.fid, '\nt simple is the simple iid t-statistic.');
  fprintf(p.fid, '\nt hetero allows for heteroscedasticity in the residual covariance matrix.');
  fprintf(p.fid, '\nt clstrd time allows for heteroscedasticity and general cross-correlation in the residual covariance matrix.');
  fprintf(p.fid, '\nt clstrd area allows for heteroscedasticity and general time-series correlation in the residual covariance matrix.');
end


function coef_stats_label = load_FM_coef_stats_label(num_autos_in_coef_se)

  coef_stats_label(1, 1) = "Coef    ";
  coef_stats_label(2, 1) = "SE simpl";
  coef_stats_label(3, 1) = "SE LHS-1";
  coef_stats_label(4, 1) = "t simpl ";
  coef_stats_label(5, 1) = "t LHS-1 ";

  for auto_indx = 1:size(num_autos_in_coef_se, 2)
    coef_stats_label(5 + auto_indx, 1) = sprintf('t %5d ', num_autos_in_coef_se(auto_indx));
  end

end


function coef_stats_label = load_Panel_coef_stats_label

  %  If x and y are not demeaned...intercept included in panel
  %     coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both

  coef_stats_label = ["Coef          "; "SE iid        "; "SE White htero"; "SE Clstrd time"; "SE Clstrd area"; "SE Clstrd both";  ...
                                        "t iid         "; "t White hetero"; "t Clstrd time "; "t Clstrd area "; "t Clstrd both "];
end 


function TS_coef_stats_label = load_TS_coef_stats_label

  TS_coef_stats_label(1, 1) = "Coef        ";
  TS_coef_stats_label(2, 1) = "SE iid      ";
  TS_coef_stats_label(3, 1) = "SE resid cor";
  TS_coef_stats_label(4, 1) = "t iid       ";
  TS_coef_stats_label(5, 1) = "t resid cor ";

end

function TS_area_coef_stats_label = load_TS_area_coef_stats_label(p)

  TS_area_coef_stats_label(1, 1) = "Coef      ";
  TS_area_coef_stats_label(2, 1) = "SE iid    ";
  TS_area_coef_stats_label(3, 1) = "SE LHS-1  ";
  TS_area_coef_stats_label(4, 1) = "SE hetero ";
  TS_area_coef_stats_label(5, 1) = "SE hac L-1";
  TS_area_coef_stats_label(6, 1) = "t iid     ";
  TS_area_coef_stats_label(7, 1) = "t LHS-1   ";
  TS_area_coef_stats_label(8, 1) = "t hetero  ";
  TS_area_coef_stats_label(9, 1) = "t hac L-1 ";

  for auto_indx = 1:size(p.num_autos_in_TS_coef_se, 2)
    TS_area_coef_stats_label(9 + auto_indx, 1) = sprintf('t %5d', p.num_autos_in_TS_coef_se(auto_indx));
  end

end

function num_coef_stats = get_num_FE_TS_coef_stats(p)   % num_coef_stats is for the CS ave of multiple TS regressions, num_area_coef_stats is for one TS regression.

  num_coef_stats  = 9;                                           % coef, SE iid, SE resid cor, t iid, t resid cor
end

function num_coef_stats = get_num_FE_CS_coef_stats(p)   % num_coef_stats is for the CS ave of multiple TS regressions, num_area_coef_stats is for one TS regression.

  num_coef_stats  = 3;                                           % coef, SE iid, t iid
end

function [num_coef_stats, num_area_coef_stats, num_ave_ts, num_area_ts] = get_num_TS_coef_stats(p)   % num_coef_stats is for the CS ave of multiple TS regressions, num_area_coef_stats is for one TS regression.

  if max(p.num_autos_in_TS_coef_se) == 0
    num_area_ts = 4;                                             % t-simple, t-LHS-1, diag(white_cov_TS)', diag(area_clstrd_cov_TS)'
  else
    num_area_ts = 4 + size(p.num_autos_in_TS_coef_se, 2);        % t-simple, t-LHS-1, diag(white_cov_TS)', diag(area_clstrd_cov_TS)', other t's adj for autocorrelation
  end

  num_area_coef_stats = 5 + num_area_ts;                         % coef, SE iid, SE LHS-1, hetero, area clustered, ts
  num_coef_stats  = 5;                                           % coef, SE iid, SE resid cor, t iid, t resid cor
  num_ave_ts      = 2;
end

function [num_coef_stats, num_addl_ts, num_autos] = get_num_FM_coef_stats(p)   % num_coef_stats = 5 + num_addl_ts = coef, SE simpl, SE LHS-1, t-simple, t-LHS-1, addl ts (e.g., 12, 24, 36)

  if max(p.num_autos_in_FM_coef_se) == 0
    num_addl_ts = 0;
  else
    num_addl_ts = size(p.num_autos_in_FM_coef_se, 2);        % other t's adj for autocorrelation...e.g. 12, 24, 36
  end

  num_coef_stats = 5 + num_addl_ts;                          % coef, SE simpl, SE LHS-1, t simple, t LHS-1  plus other t's adj for autocorrelation...e.g. 12, 24, 36

  num_autos = max([1, p.LHS_per - 1, p.num_autos_FM_coefs, p.num_autos_in_FM_coef_se]);

end


function [num_coef_stats, num_ts] = get_num_mthly_dummies_coef_stats
  num_ts         = 2;              % t-simple, t-y_per - 1
  num_coef_stats = 3 + num_ts;     % coef, SE simple, SE LHS_per - 1, ts
end

function num_coef_stats = get_num_area_dummies_coef_stats
    num_coef_stats = 5;     % coef, SE simple, SE resid cor, t simple, t resid cor
end

function coef_stats_label = load_dummies_coef_stats_label(num_autos_in_coef_se)

  coef_stats_label(1, 1) = "Coef    ";
  coef_stats_label(2, 1) = "SE simpl";
  coef_stats_label(3, 1) = "SE LHS-1";
  coef_stats_label(4, 1) = "t simpl ";
  coef_stats_label(5, 1) = "t LHS-1 ";
end

function num_coef_stats = get_num_Panel_coef_stats

  %     coef_stats = [P'; SE_P; t_stats];   % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, clstrd area, clstrd both;  t-stats: iid, hetero, clstrd time, area, both
  %
  %
  num_coef_stats = 11;                              % coef, SEs: iid, hetero, clstrd time, area, both               t-stats: iid, hetero, clstrd time, area, both
end

function x_adj = Panel_demean_x(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x_raw, style_passed, c, d, p)

  if style_passed == c.demean_both_full
    x_adj = demean_x_both_full(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x_raw, c, d, p);
  
  elseif x_in_y(style_passed, [c.demean_CS_w_int, c.demean_CS_wo_int, c.demean_both_CS_1st_w, c.demean_both_CS_1st_wo])
    x_adj  = demean_x_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x_raw, c, d, p);
 
    if x_in_y(style_passed, [c.demean_both_CS_1st_w, c.demean_both_CS_1st_wo])
      x_adj  = demean_x_each_area(frst_form_mth_area, last_form_mth_area, x_adj, c, d, p);
    end

  elseif x_in_y(style_passed, [c.demean_TS_w_int, c.demean_TS_wo_int, c.demean_both_TS_1st_w, c.demean_both_TS_1st_wo])
    x_adj  = demean_x_each_area(frst_form_mth_area, last_form_mth_area, x_raw, c, d, p);
 
    if x_in_y(style_passed, [c.demean_both_TS_1st_w, c.demean_both_TS_1st_wo])
      x_adj  = demean_x_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x_adj, c, d, p);
    end
  
  else    % x_adj is x_raw and LHS_adj is LHS_raw if not demeaning
    x_adj  = x_raw;
  end
end

function y_adj = Panel_demean_y(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y_raw, style_passed, c, d, p)

  if style_passed == c.demean_both_full
    y_adj = demean_y_both_full(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y_raw, c, d, p);

  elseif x_in_y(style_passed, [c.demean_CS_w_int, c.demean_CS_wo_int, c.demean_both_CS_1st_w, c.demean_both_CS_1st_wo])
    y_adj = demean_y_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y_raw, c, d, p);

    if x_in_y(style_passed, [c.demean_both_CS_1st_w, c.demean_both_CS_1st_wo])
      y_adj = demean_y_each_area(frst_form_mth_area, last_form_mth_area, LHS_per, y_adj, c, d, p);
    end

  elseif x_in_y(style_passed, [c.demean_TS_w_int, c.demean_TS_wo_int, c.demean_both_TS_1st_w, c.demean_both_TS_1st_wo])
    y_adj = demean_y_each_area(frst_form_mth_area, last_form_mth_area, LHS_per, y_raw, c, d, p);

    if x_in_y(style_passed, [c.demean_both_TS_1st_w, c.demean_both_TS_1st_wo])
      y_adj = demean_y_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y_adj, c, d, p);
    end
  
  else    % PR_adj is PR_raw and y_adj is y_raw if not demeaning
    y_adj = y_raw;
  end
end

function x_dmnd = demean_x_each_area(frst_form_mth_area, last_form_mth_area, x, c, d, p)

  x_dmnd = NaN(size(x));
  
  % subtract TS means for each area
  %
  for area = 1:c.num_areas_both
    
    frst_x = frst_form_mth_area(area);  
    last_x = last_form_mth_area(area);  

    if ismatrix(x) 
      x_dmnd(frst_x:last_x, area) = x(frst_x:last_x, area) - mean(x(frst_x:last_x, area));
    else
      x_dmnd(frst_x:last_x, :, area) = x(frst_x:last_x, :, area) - mean(x(frst_x:last_x, :, area));
    end
  end

end

function y_dmnd = demean_y_each_area(frst_form_mth_area, last_form_mth_area, LHS_per, y, c, d, p)

  y_dmnd = NaN(size(y));

  % subtract TS means for each area
  %
  % y is observed LHS_per months after formation mth
  %
  for area = 1:c.num_areas_both
    
    frst_y = frst_form_mth_area(area) + LHS_per; 
    last_y = last_form_mth_area(area) + LHS_per; 
    
    y_dmnd(frst_y:last_y, area) = y(frst_y:last_y, area) - mean(y(frst_y:last_y, area), 'omitnan');
  end

end

function y_dmnd = demean_y_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y, c, d, p)

  % load matrices with just good months for each area...values outside frst and last months will be NaN
  %
  y_NaN  = NaN(size(y));
  y_dmnd = NaN(size(y));

  for area = 1:c.num_areas_both
    frst_y = max(frst_form_mth, frst_form_mth_area(area)) + LHS_per;
    last_y = min(last_form_mth, last_form_mth_area(area)) + LHS_per;

    y_NaN(frst_y:last_y, area) = y(frst_y:last_y, area);
  end
  
  for mth_y = frst_form_mth + LHS_per:last_form_mth +LHS_per
    y_dmnd(mth_y, :) = y_NaN(mth_y, :) - mean(y_NaN(mth_y, :), 'omitnan');
  end

end

function x_dmnd = demean_x_each_month(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x, c, d, p)

  % load matrices with just good months for each area...values outside frst and last months will be NaN
  %
  x_NaN  = NaN(size(x));
  x_dmnd = NaN(size(x));

  for area = 1:c.num_areas_both
    frst_x = max(frst_form_mth, frst_form_mth_area(area)); 
    last_x = min(last_form_mth, last_form_mth_area(area)); 

    if ismatrix(x) 
      x_NaN(frst_x:last_x, area) = x(frst_x:last_x, area);
    else
      x_NaN(frst_x:last_x, :, area) = x(frst_x:last_x, :, area);
    end
  end
  
  for mth_x = frst_form_mth:last_form_mth
    if ismatrix(x) 
      x_dmnd(mth_x, :) = x_NaN(mth_x, :) - mean(x_NaN(mth_x, :), 'omitnan');
    else
      x_dmnd(mth_x, :, :) = x_NaN(mth_x, :, :) - mean(x_NaN(mth_x, :, :), 'omitnan');
    end
  end

end


function x_dmnd = demean_x_both_full(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, x, c, d, p)

  % load matrices with just good months for each area...values outside frst and last months will be NaN
  %
  x_NaN  = NaN(size(x));
  x_dmnd = NaN(size(x));

  for area = 1:c.num_areas_both
    frst_x = max(frst_form_mth, frst_form_mth_area(area)); 
    last_x = min(last_form_mth, last_form_mth_area(area)); 

    if ismatrix(x) 
      x_NaN(frst_x:last_x, area) = x(frst_x:last_x, area);
    else
      x_NaN(frst_x:last_x, :, area) = x(frst_x:last_x, :, area);
    end
  end
  
  % subtract cross-section mean for each month
  %
  for mth_x = frst_form_mth:last_form_mth
    if ismatrix(x) 
      x_dmnd(mth_x, :) = x_NaN(mth_x, :) - mean(x_NaN(mth_x, :), 'omitnan');
    else
      x_dmnd(mth_x, :, :) = x_NaN(mth_x, :, :) - mean(x_NaN(mth_x, :, :), 'omitnan');
    end
  end


  % subtract TS means for each area
  %
  for area = 1:c.num_areas_both
    
    frst_x = frst_form_mth_area(area);  
    last_x = last_form_mth_area(area);  

    if ismatrix(x) 
      x_dmnd(frst_x:last_x, area) = x_dmnd(frst_x:last_x, area) - mean(x(frst_x:last_x, area), 'omitnan');
    else
      x_dmnd(frst_x:last_x, :, area) = x_dmnd(frst_x:last_x, :, area) - mean(x(frst_x:last_x, :, area), 'omitnan');
    end
  end


  % subtract grand mean, so average x_dmnd is zero
  %
  if p.subtract_grnd_mn_both_full
    if ismatrix(x) 
      x_dmnd = x_dmnd - mean(x_dmnd, 'all', 'omitnan');
    else
      for x_var = 1:size(x, 2)
        x_dmnd(:, x_var, :) = x_dmnd(:, x_var, :) - mean(x_dmnd(: , x_var, :), 'all', 'omitnan');
      end
    end
  end
end


function x_dmnd = demean_x_both_smpl(x, c, d, p)

  % load matrices with just good years for each area...values outside frst and last years will be NaN
  %
  num_yrs  = size(x, 1);
  num_MSAs = size(x, 2);
  num_vars = size(x, 3);
  
  x_dmnd = NaN(size(x));

  for var = 1:num_vars
  
    % subtract cross-section mean for each year
    %
    for yr = 1:num_yrs
      x_dmnd(yr, :, var) = x(yr, :, var) - mean(x(yr, :, var), 'omitnan');
    end
    
    % subtract TS means for each area
    %
    for MSA = 1:num_MSAs
      x_dmnd(:, MSA, var) = x_dmnd(:, MSA, var) - mean(x(:, MSA, var), 'omitnan');
    end
    
    % subtract grand mean, so average var_dmnd is zero
    %
    if p.subtract_grnd_mn_both_full
      x_dmnd(:, :, var) = x_dmnd(:, :, var) - mean(x_dmnd(:, :, var), 'all', 'omitnan');
    end
        
  end
end


function y_dmnd = demean_y_both_full(frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, LHS_per, y, c, d, p)

  % load matrices with just good months for each area...values outside frst and last months will be NaN
  %
  y_NaN  = NaN(size(y));
  y_dmnd = NaN(size(y));

  for area = 1:c.num_areas_both
    frst_y = max(frst_form_mth, frst_form_mth_area(area)) + LHS_per; 
    last_y = min(last_form_mth, last_form_mth_area(area)) + LHS_per; 

    y_NaN(frst_y:last_y, area) = y(frst_y:last_y, area);
  end
  
  % subtract cross-section mean for each month
  %
  for mth_y = frst_form_mth + LHS_per:last_form_mth + LHS_per
    y_dmnd(mth_y, :) = y_NaN(mth_y, :) - mean(y_NaN(mth_y, :), 'omitnan');
  end


  % subtract TS means for each area
  %
  for area = 1:c.num_areas_both
    
    frst_y = frst_form_mth_area(area) + LHS_per;  
    last_y = last_form_mth_area(area) + LHS_per;  

    y_dmnd(frst_y:last_y, area) = y_dmnd(frst_y:last_y, area) - mean(y(frst_y:last_y, area), 'omitnan');
  end


  % subtract grand mean, so average y_dmnd is zero
  %
  if p.subtract_grnd_mn_both_full
    y_dmnd = y_dmnd - mean(y_dmnd, 'all', 'omitnan');
  end
  

end

function y_dmnd = demean_y_both_smpl(y, c, d, p)

  % load matrices with just good years for each area...values outside frst and last years will be NaN
  %
  num_mths = size(y, 1);
  num_MSAs = size(y, 2);
  
  y_dmnd = NaN(size(y));

  % subtract cross-section mean for each year
  %
  for mth = 1:num_mths
    y_dmnd(mth, :) = y(mth, :) - mean(y(mth, :), 'omitnan');
  end
    
  % subtract TS means for each area
  %
  for MSA = 1:num_MSAs
    y_dmnd(:, MSA) = y_dmnd(:, MSA) - mean(y(:, MSA), 'omitnan');
  end
    
  % subtract grand mean, so average var_dmnd is zero
  %
  if p.subtract_grnd_mn_both_full
    y_dmnd = y_dmnd - mean(y_dmnd, 'all', 'omitnan');
  end
        
end


function y_dmnd = demean_y_TS_smpl(y, c, d, p)

  % load matrices with just good years for each area...values outside frst and last years will be NaN
  %
  num_MSAs = size(y, 2);
  
  y_dmnd = NaN(size(y));

  % subtract TS means for each area
  %
  for MSA = 1:num_MSAs
    y_dmnd(:, MSA) = y(:, MSA) - mean(y(:, MSA), 'omitnan');
  end
end


function y_dmnd = demean_y_CS_smpl(y, c, d, p)

  % load matrices with just good years for each area...values outside frst and last years will be NaN
  %
  num_mths  = size(y, 1);

  y_dmnd = NaN(size(y));

  % subtract cross-section mean for each year
  %
  for mth = 1:num_mths
    y_dmnd(mth, :) = y(mth, :) - mean(y(mth, :), 'omitnan');
  end
end




function t_or_f = x_in_y(x, y)

  if any(y == x)
    t_or_f = true;
  else
    t_or_f = false;
  end
end


function coef_stats_out = add_sum_stats_for_implied_intcept_Panel(y_raw, x_raw, nt, source, LHS_per, P, style_passed, coef_stats_in, c, d, p)


  % add summary stats for implied intercept if not included in x
  %
  % compute monthly implied intercepts...include if demean_CS_wo_int or demean_both wo int and TS first...so CS is second
  %
  if x_in_y(style_passed, [c.demean_CS_wo_int, c.demean_both_TS_1st_wo])
  
    num_mths  = size(nt, 1);

    % Add intercept implied by CS average of x and y each month
    %
    intcept = NaN(num_mths, 1);   % num_mths x 1
    
    obs_end = 0;
    for t = 1:num_mths
      obs_beg = obs_end + 1;
      obs_end = obs_end + nt(t);
      intcept(t) = mean(y_raw(obs_beg:obs_end)) - mean(x_raw(obs_beg:obs_end, :)) * P;
    end

    [coef_stats_shrt, ~] = summarize_FM_coefs(intcept, LHS_per, c, d, p);    % (3 + 2 + num_addl_ts) x num_coefs = (ave coef, se_smpl, se_LHS_1, addl_ts... num_autos x num_coefs  ~ replaces autos

    coef_stats_out         = [NaN(size(coef_stats_in, 1), 1), coef_stats_in];
    coef_stats_out(1:5, 1) = coef_stats_shrt(1:5, 1);  % ave int; SE ave int; auto-adjusted SE_ave_int; t-simple; t LHS per-1

    
  elseif x_in_y(style_passed, [c.demean_TS_wo_int, c.demean_both_CS_1st_wo])

    % Add intercept implied by TS average for each area
    %
    intcept = NaN(c.num_areas_both, 1);   % c.num_areas_both x 1
    
    for area = 1:c.num_areas_both
      x_raw_area = x_raw(source == area, :);
      y_raw_area = y_raw(source == area);
      intcept(area) = mean(y_raw_area) - mean(x_raw_area) * P;
    end
    
    coef_stats_out               = [NaN(size(coef_stats_in, 1), 1), coef_stats_in];
    coef_stats_out([1, 2, 4], 1) = summarize_TS_area_coefs_simple(intcept, c, d, p);   % out is (ave coef, se_smpl, se_LHS_1, t_smpl, t_LHS_1), summarize_TS_area_coefs_simple produces ave coef, se_smpl, t_smpl

  elseif style_passed == c.demean_both_full
    coef_stats_out               = [NaN(size(coef_stats_in, 1), 1), coef_stats_in];
  end

end


function frst_slope = determine_frst_slope(style_passed, c, d, p)

  if x_in_y(style_passed, [c.not_demeaned, c.demean_TS_w_int, c.demean_CS_w_int, c.demean_both_TS_1st_w, c.demean_both_CS_1st_w])
    frst_slope = 2;
  else
    frst_slope = 1;
  end
end


function [ave_SD_x_adj, SD_SD_x_adj] = summarize_CS_SD_of_x_adj(x_adj, nt, style_passed, c, d, p)

  num_mths   = size(nt,    1);
  num_coefs  = size(x_adj, 2);
  frst_slope = determine_frst_slope(style_passed, c, d, p);
  num_slopes = num_coefs - frst_slope + 1;

  SD_x_adj  = zeros(num_slopes, num_mths);  % This include column of 1 if regression includes intercept...The first column SD_x_adj is dropped later

  obs_end   = 0;
  for t = 1:num_mths
    obs_beg = obs_end + 1;           
    obs_end = obs_end + nt(t);             % nt is number of areas in month t

    x_adj_t = x_adj(obs_beg:obs_end, :);   % nt x num_coefs
        
    SD_x_adj(:, t) = std(x_adj_t(:, frst_slope:num_coefs))';  % drop intercept if there is one
  end 

  % Std_x_adj does not include column for intercept
  %
  ave_SD_x_adj = mean(SD_x_adj, 2);       % num_slopes x 1                       
  SD_SD_x_adj  =  std(SD_x_adj, 0, 2);    % num_slopes x 1       0 is so it uses normal (equal) weighting function and divides by n-1

end


% ******************************************************************************
%                                                                              *
%                         Functions for All Regressions                      *
%                                                                              *
% ******************************************************************************

function [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(FM_TS_or_Panel, style_passed, c, d, p)

  num_RHS_pers       = size(p.RHS_per, 2);
  max_RHS_per        =  max(p.RHS_per);
  num_LHS_pers       = size(p.LHS_per, 2);
  num_LHS_RHS_combos = size(p.LHS_RHS_vars, 2);
  num_regrs_base     = 1 + num_RHS_pers;      % PR, PR & dP(t-K,t)
  num_regrs          = num_regrs_base;
    
  max_num_slopes     = 2;
  max_num_coefs      = 3;
  
  frst_slope         = determine_frst_slope(style_passed, c, d, p);   % frst_slope = 2 if there is an intercept
  
  if FM_TS_or_Panel == c.FM
    num_autos      = p.num_autos_FM_coefs;
  elseif FM_TS_or_Panel == c.TS
    num_autos      = max([p.num_autos_resids, p.num_autos_in_TS_coef_se]);
  elseif FM_TS_or_Panel == c.FE_TS
    num_autos      = max([p.num_autos_FE_resids, p.LHS_per - 1]);
  elseif FM_TS_or_Panel == c.Panel
    num_autos      = 1;                  % not really used
  elseif FM_TS_or_Panel == c.FE_CS
    num_autos      = 1;                  % not really used
  else
    disp('FM_TS_or_Panel must be c.FM, c.TS, c.Panel, c.FE_TS, or c.FE_CS in setup_wrt_and_run_regr_constants')
    pause
  end

end


function TS_reg_results = TS_FE_regr(x, y, y_per, c, d, p)

  % y is num_obs x 1
  % x is num_obs x num_coefs
  %
  % coef_stats is 7 x num_coefs:  coef, iid SE, hetero SE, area clstrd SE, t iid, t hetero, t area clustrd
  % glbl_r2_ttl_raw is scalar
  num_obs   = size(x, 1);
  num_coefs = size(x, 2);
  
  
  % Compute coefficients
  %
  coef = (x' * x) \ x' * y;
  %Kx1    KxT TxK   KxT Tx1
  
  coef = coef';
  %1xK
  

  % The CS average for each month may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y = var(y);


  % Calculate residuals and adjust for overfitting if demeaned
  %
  resids = y - x * coef';
  %       Tx1 TxK   Kx1
           

  % Calculate residual variances and r2
  %
  var_e = sum(resids.^2, 1) / (num_obs - num_coefs);
  r2    = 1.0 - (var_e ./ var_y);          % 1 x 1

  % Compute SEs of coefs assuming iid resids
  %
  xx     = x' * x;                                     % num_coefs x num_coefs
  inv_xx = inv(xx);

  iid_cov_TS = var_e * inv_xx;
  %   KxK       1x1    K x K


  % sum autos up to y_per - 1 
  %
  if y_per == 1
    auto_adj = 1;
  else
    num_lags   = y_per - 1;
    autos_shrt = autocorr(resids, 'NumLags', num_lags);
    autos      = autos_shrt(2:num_lags + 1, 1);    % num_autos x 1

    auto_adj   = cmpt_auto_adj(num_lags, num_obs, autos(1:num_lags));
  end

  
  % Heteroscedastic consistent SEs
  %
%  disp("Hetero")
%  z_tmp = diag(resids_adj.^2);
%  size(z_tmp)
%  z_tmp(1:10, 1:10)
  
  x_omega_x_white = x' * diag(resids.^2) * x;          % diag makes a square matrix with resid_adj^2 on the diagonal

  white_cov_TS = inv_xx * x_omega_x_white * inv_xx;    % num_coefs x num_coefs  
  

  % Generalized HH
  %
  sig2 = resids * resids';
  
  for row = 1:num_obs
    for col = 1:row - y_per
      sig2(row, col) = 0;
      sig2(col, row) = 0;
    end
  end
  
  genrlzd_HH_x_sig2_x = x' * sig2 * x;
  %                    KxT  T x T  TxK
      
  genrlzd_HH_cov_TS = inv_xx * genrlzd_HH_x_sig2_x * inv_xx;
  
      SE_TS = NaN(4, num_coefs);
  t_stat_TS = NaN(4, num_coefs);
  
  SE_TS(1:4, :) = [diag(iid_cov_TS)'; auto_adj * diag(iid_cov_TS)'; diag(white_cov_TS)'; diag(genrlzd_HH_cov_TS)'];
  % 4xK
  
  SE_TS = sqrt(SE_TS);

  for indx = 1:size(SE_TS, 1)
    t_stat_TS(indx, :) = coef ./ SE_TS(indx, :);
  end
  
  TS_reg_results = NaN(10, num_coefs);
  TS_reg_results(1:9, 1:num_coefs) = [coef; SE_TS; t_stat_TS];         % coefs, se_coefs, t-stats  1 + 4 + 4
  TS_reg_results(10,  1)           = r2;
  
end


%-----------------------------------------------------------------------------------------
%                                                                                        *
%                                          FM Regressions                                *
%                                                                                        *
% ****************************************************************************************
%
function FM = run_FM_dR_on_PR_and_one_dP(style_passed, c, d, p)

  % Runs FM regressions of dR on PR and one dP. Called by Main.
  % "One" means for each combination of gap and dR, we run just PR, then PR and (one by one) something like the 1, 3, 6, and 12 mth change in dP to explain dR.
  % Calls run_set_of_FM_dR_on_PR_and_one_dP for each combiation of gap, LHS_per, and RHS_per. The function computes regrs for all months. 
  %   The subfunction run_set_of_FM_dR_on_PR_and_one_dP calls summarize_FM_coefs and returns the results.
  %
  % This function runs FM regressions of dR or dP on two or three variables -- such as dR on intercept, P/R, and dP, or dR on P/R and dP, or dP on P/R and dP
  % THe LHS variable might be the change in R from formation month t, when the RHS variables are observed to t+1, 3, 6, or 12.
  % The RHS variable might be the change in P from t-1 to t, t-3 to t, t-6 to t, or t-12 to t.
  %
  % c = area
  % t = month
  % e2ct = squared residual for area c in month t
  % Nt = number of areas in month t
  % Nc = number of months for area c
  % K  = number of regression coefficients
  %
  % style_passed
  %
  % c.not_demeaned                   FM       SE2t = {sum c e2ct} / (Nt - K)                        SE2 = sum t {Nt SE2t / N}
  %								  
  % c.demean_TS_w_int or _wo_int     FM       SE2t = {sum c [Nc/(Nc - 1)]e2ct} / (Nt - K)           SE2 = sum t {Nt SE2t / N}
  %
  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.FM, style_passed, c, d, p);
  
  [num_coef_stats, num_addl_ts, num_autos] = get_num_FM_coef_stats(p);                  % num_coef_stats = 5 + num_addl_ts = coef, SE simpl, SE LHS-1, t-simple, t-LHS-1, addl ts (e.g., 12, 24, 36)

  coef_stats_label          =  load_FM_coef_stats_label(p.num_autos_in_FM_coef_se);                                
										 
  FM.coef_stats             =  NaN(num_coef_stats, max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % simple stats + num_addl_ts
  FM.glbl_r2_ttl_raw        =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % 1 - variance of bias adjusted residuals relative to variance of raw y
  FM.glbl_r2_dmn_raw        =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y
  FM.glbl_r2_slps_raw_y          =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by x = R2 total - R2 demean
  FM.glbl_r2_slps_dmn_y          =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y
  FM.glbl_var_e             =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.glbl_var_raw_y         =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.glbl_var_BA_adj_y      =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  FM.FE_r2                  =                                          zeros(num_LHS_pers, num_LHS_RHS_combos);
  FM.ave_r2                 =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.ave_var_e              =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.autos                  =  NaN(num_autos,      max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.ave_SD_x_adj           =               zeros(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.SD_SD_x_adj            =               zeros(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.cov_PR_LHS             =                          NaN(num_RHS_pers + 1, num_LHS_pers, num_LHS_RHS_combos);
  FM.cor_PR_dP              =                              NaN(num_RHS_pers, num_LHS_pers, num_LHS_RHS_combos);
  FM.ave_resid              =               NaN(c.num_areas_both, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.t_resid                =               NaN(c.num_areas_both, num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  FM.glbl_cov_pnl_wght_coef =                  NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.glbl_cor_pnl_wght_coef =                  NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FM.glbl_cor_var_x_coef    =                  NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  % Cycle over gaps etc
  %
  for LHS_RHS_var_indx = 1:num_LHS_RHS_combos
    LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
    RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);
      
    for LHS_per_indx = 1:num_LHS_pers
      LHS_per = p.LHS_per(LHS_per_indx);

      FM_tmp = run_set_of_FM_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p);

               FM.coef_stats(:, :, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.coef_stats;
                FM.glbl_r2_ttl_raw(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_r2_ttl_raw;
                FM.glbl_r2_dmn_raw(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_r2_dmn_raw;
             FM.glbl_r2_slps_raw_y(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_r2_slps_raw_y;
             FM.glbl_r2_slps_dmn_y(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_r2_slps_dmn_y;

                     FM.glbl_var_e(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_var_e;
                 FM.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_var_raw_y;
              FM.glbl_var_BA_adj_y(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_var_BA_adj_y;

      FM.glbl_cov_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_cov_pnl_wght_coef;
      FM.glbl_cor_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_cor_pnl_wght_coef;
         FM.glbl_cor_var_x_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.glbl_cor_var_x_coef;

                         FM.ave_r2(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.ave_r2;
                      FM.ave_var_e(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.ave_var_e;
                FM.ave_SD_x_adj(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.ave_SD_x_adj;     % This is the average of the CS standard deviations of the demeaned Xs
                 FM.SD_SD_x_adj(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.SD_SD_x_adj;      % This is the std dev of the CS standard deviations of the demeaned Xs
                    FM.autos(:, :, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.autos;
                     FM.cov_PR_LHS(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.cov_PR_LHS;
                      FM.cor_PR_dP(:, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.cor_PR_dP;
                   FM.ave_resid(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.ave_resid;
                     FM.t_resid(:, :, LHS_per_indx, LHS_RHS_var_indx) = FM_tmp.t_resid;
    end
  end
end

% OK
%
function FM_tmp = run_set_of_FM_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p)

  % Set here means all regressions for combination of gap and LHS_per, we run just PR, then PR and (one by one) the 1, 3, 6, and 12 mth change in dP to explain dR for all months.
  % Computes coefs and r2s directly in the function.
  % Uses summarize_FM_coefs to compute average coefs etc.
  % Has the option to write the monthly coefficients that contribute the most to each average coefficient's t-stat.

  %  coef_stats  = zeros(num_coef_stats, max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  %  autos       = zeros(num_autos,      max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  %  glbl_r2     =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos); 
  %  glbl_var_e  =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos); 
  %  ave_r2      =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos); 
  %  ave_var_e   =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos); 

  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.FM, style_passed, c, d, p);


  [num_coef_stats, num_addl_ts, num_autos] = get_num_FM_coef_stats(p);                  % num_coef_stats = 5 + num_addl_ts = coef, SE simpl, SE LHS-1, t-simple, t-LHS-1, addl ts (e.g., 12, 24, 36)

  % Global results
  %
  coef_stats             = NaN(num_coef_stats,  max_num_coefs,  num_regrs);
  autos                  = NaN(num_autos,       max_num_slopes, num_regrs);
				
  ave_SD_x_adj           =                  NaN(max_num_slopes, num_regrs);
  SD_SD_x_adj            =                  NaN(max_num_slopes, num_regrs);
  glbl_r2_ttl_raw        =                                  NaN(num_regrs, 1);    % 1 - variance of bias adjusted residuals relative to variance of raw y = 1 - var_e/var_raw_y)
  glbl_r2_dmn_raw        =                                  NaN(num_regrs, 1);    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
  glbl_r2_slps_raw_y     =                                  NaN(num_regrs, 1);    % fraction of raw variance explained by x = R2 total - R2 demean = 
  glbl_r2_slps_dmn_y     =                                  NaN(num_regrs, 1);    % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y

  glbl_var_e             =                                zeros(num_regrs, 1);
  glbl_var_raw_y         =                                zeros(num_regrs, 1);
  glbl_var_BA_adj_y      =                                zeros(num_regrs, 1);
  
  glbl_cov_pnl_wght_coef = NaN(max_num_coefs, num_regrs);
  glbl_cor_pnl_wght_coef = NaN(max_num_coefs, num_regrs);
  glbl_cor_var_x_coef    = NaN(max_num_coefs, num_regrs);
  
  ave_r2                 = NaN(num_regrs, 1);
  ave_var_e              = NaN(num_regrs, 1);
  ave_cov_PR_LHS         = NaN(num_RHS_pers + 1, 1);
  ave_cor_PR_dP          = NaN(num_RHS_pers, 1);
  ave_resid              = NaN(c.num_areas_both, num_regrs);
  t_resid                = NaN(c.num_areas_both, num_regrs);
  
  % Working arrays
  %
  coef_mth      = NaN(c.last_mth, max_num_coefs);           % If no intercept, will use implied intercept
  r2_mth        = NaN(c.last_mth, 1);
  var_e_mth     = NaN(c.last_mth, 1);
  SD_x_adj_mth  = NaN(c.last_mth, max_num_slopes);

  y_raw_ttl     = NaN(c.max_TS_CS_obs, 1);
  y_BA_adj_ttl  = NaN(c.max_TS_CS_obs, 1);

  y_raw_mth     = NaN(c.num_areas_both, 1);
  y_adj_mth     = NaN(c.num_areas_both, 1);
  y_BA_adj_mth  = NaN(c.num_areas_both, 1);
  x_raw_mth     = NaN(c.num_areas_both, max_num_coefs);
  x_adj_mth     = NaN(c.num_areas_both, max_num_coefs);
  resids_FE_mth = NaN(c.num_areas_both, 1);
   srce_mth     = NaN(c.num_areas_both, 1);

  sum_x2        = NaN(c.last_mth, max_num_coefs);
  var_x         = NaN(c.last_mth, max_num_coefs);
  pnl_wght      = NaN(c.last_mth, max_num_coefs);


  if style_passed == c.not_demeaned
    x_adj_mth(:, 1) = 1;
    x_raw_mth(:, 1) = 1;
  end

  % determine and possibly write frst and last form mth
  %
  % frst_form_mth = NaN(num_regrs_base, 1);  
  % last_form_mth = NaN;                       
  % frst_form_mth_area = NaN(c.num_areas_both, num_regrs_base);    % 1 + num_RHS_pers...1 is for simple P/R
  % last_form_mth_area = NaN(c.num_areas_both, 1);            % 1 + num_RHS_pers
  %
  % nt = zeros(c.num_mths_ttl,    num_RHS_pers + 1)            % number of good areas each month
  % nc = zeros(c.num_areas_both, num_RHS_pers + 1)            % number of good months for each area
  %
  if LHS_per == 1 && ~p.wrt_limited_output
    fprintf(p.fid, '\n\nFM  LHS_per = %d  LHS_var = %d  RHS_var = %d', LHS_per, LHS_var, RHS_var);
    [frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, nt, nc] = determine_and_wrt_frst_and_last_form_mth(LHS_per, LHS_var, RHS_var, c.wrt_frst_last_mth, c, d, p);
  else
    [frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, nt, nc] = determine_and_wrt_frst_and_last_form_mth(LHS_per, LHS_var, RHS_var, c.dont_wrt_frst_last_mth, c, d, p);
  end


  % Select LHS_chg data based on LHS_var and LHS_per...0 means for all areas
  %
  LHS_raw = load_mthly_qtrly_semi_or_annual_chg(LHS_var, LHS_per, 0, c, d, p);  % The two omitted results are for subtracting the overall mean for the month (in the TS regression)
  

  % Load vector of PR or RP data
  %
  PR_raw = load_full_PR_vector_regr(0, c, d, p);

  if style_passed == c.not_demeaned
    LHS_adj    = LHS_raw;
    PR_adj     = PR_raw;
    resids_FE  = LHS_raw;
  end
  
  % Run FM regression with only PR if regr = 1 or PR and dP if regr > 1
  %
  for regr = 1:num_regrs

    resid_cnt    = zeros(c.num_areas_both, 1);
    resid_mth    = NaN(last_form_mth, c.num_areas_both);

    % Setup RHS variable
    %
    if regr == 1
      last_slope = frst_slope;
      num_slopes = 1;
      regr_base  = regr;
      cov_PR_LHS = NaN(c.last_mth, 1);
      
    else
      last_slope = frst_slope + 1;
      num_slopes = 2;
      regr_base  = regr;
      cov_PR_LHS = NaN(c.last_mth, 1);
      cor_PR_dP  = NaN(c.last_mth, 1);
    end
    

    if regr ~= 1
      RHS_per      = p.RHS_per(regr_base - 1);
      RHS_raw_lag  = load_mthly_qtrly_semi_or_annual_chg(RHS_var, RHS_per, 0, c, d, p);  % This loads lagged RHS_var for all areas
      RHS_raw_lead = load_mthly_qtrly_semi_or_annual_chg(RHS_var, LHS_per, 0, c, d, p);  % This loads future RHS_var for all areas
      
      if style_passed == c.not_demeaned
        RHS_adj_lag = RHS_raw_lag;
      elseif style_passed == c.demean_CS_wo_int
        RHS_adj_lag = demean_x_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, RHS_raw_lag, c, d, p);
      end
    end


    if style_passed == c.demean_CS_wo_int
      LHS_adj = demean_y_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, LHS_raw, c, d, p);
      PR_adj  = demean_x_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, PR_raw, c, d, p);

      % adjust resids_FE for CS demeaning bias
      %
      if regr == 1
        resids_FE = LHS_adj;
        for mth = frst_form_mth(regr_base):last_form_mth
          resids_FE(mth, :) = resids_FE(mth, :) * sqrt(nt(mth, regr) / (nt(mth, regr) - 1));
        end
      end

    
    elseif style_passed ~= c.not_demeaned
      LHS_adj = demean_y_each_area(frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, LHS_raw, c, d, p);
      PR_adj  = demean_x_each_area(frst_form_mth_area(:, regr_base), last_form_mth_area, PR_raw, c, d, p);
      
      % adjust resids_FE for TS demeaning bias
      %
      if regr == 1
        resids_FE = LHS_adj;
        for area = 1:c.num_areas_both
          resids_FE(:, area) = resids_FE(:, area) * sqrt(nc(area, regr) / (nc(area, regr) - 1));
        end
      end
    end
    
        
    % Compute FM regressions with only PR if regr = 1, PR and dP-K if 1 < regr
    %
    SSE         = 0;
    SSE_FE      = 0;
    num_obs_ttl = 0;

    for mth = frst_form_mth(regr_base):last_form_mth
      num_obs = 0;
      for area = 1:c.num_areas_both
        if (0 < min(d.frst_ob(area, :))) && (frst_form_mth_area(area, regr_base) <= mth) && (mth <= last_form_mth_area(area))
          num_obs                          = num_obs + 1;
          y_raw_mth(num_obs)               = LHS_raw(mth + LHS_per, area);
          y_adj_mth(num_obs)               = LHS_adj(mth + LHS_per, area);

          if style_passed == c.not_demeaned
            y_BA_adj_mth(num_obs)          = LHS_adj(mth + LHS_per, area);
            
          else   % style_passed = c.demean_CS_wo_int
            y_BA_adj_mth(num_obs)          = LHS_adj(mth + LHS_per, area) * sqrt(nt(mth, regr_base) / (nt(mth, regr_base) - 1));
          end
                      
          x_raw_mth(num_obs, frst_slope) = PR_raw(mth, area);
          x_adj_mth(num_obs, frst_slope) = PR_adj(mth, area);
		
          if 1 < regr
            x_raw_mth(num_obs, frst_slope + 1) = RHS_raw_lag(mth, area);
            x_adj_mth(num_obs, frst_slope + 1) = RHS_adj_lag(mth, area);
          end

          srce_mth(num_obs) = area;
	
          % Store FE resids for month...These are bias-adjusted if TS demeaned
          %
          if regr == 1
            resids_FE_mth(num_obs) = resids_FE(mth + LHS_per, area);
          end          
        end
      end
      
      % Subtract monthly mean, adjust for bias, and cumulate SSE_FE
      % 
      if regr == 1
        resids_FE_mth(1:num_obs) = (resids_FE_mth(1:num_obs) - mean(resids_FE_mth(1:num_obs)));
        SSE_FE = SSE_FE + (num_obs / (num_obs - 1)) * sum(resids_FE_mth(1:num_obs).^2, 1);
      end          
      
      
      cov_mat         = cov(x_adj_mth(1:num_obs, frst_slope), y_adj_mth(1:num_obs));
      cov_PR_LHS(mth) = cov_mat(1, 2);
        
      if 1 < regr
        cor_PR_dP(mth)  = corr(x_adj_mth(1:num_obs, frst_slope), x_adj_mth(1:num_obs, last_slope));
      end
	
	
      SD_x_adj_mth(mth, 1:num_slopes) = std(x_adj_mth(1:num_obs, frst_slope:last_slope))';
      
      [coef_mth(mth, 1:last_slope), r2_mth(mth), var_e_mth(mth), resids_adj, y_hat_mth] = cmpt_FM_coefs_and_r2(x_adj_mth(1:num_obs, 1:last_slope), y_raw_mth(1:num_obs), y_adj_mth(1:num_obs), srce_mth(1:num_obs), nc(:, regr_base), style_passed, c, d, p);
      
      for obs = 1:num_obs
        resid_cnt(srce_mth(obs), 1)                        = resid_cnt(srce_mth(obs), 1) + 1;
        resid_mth(resid_cnt(srce_mth(obs)), srce_mth(obs)) = resids_adj(obs);
      end

      
      % Cumulate SSE and store y_raw for all observations
      %
      if style_passed == c.demean_CS_wo_int
        SSE = SSE + (num_obs / (num_obs - last_slope - 1)) * sum(resids_adj.^2, 1);
      else
        SSE = SSE + (num_obs / (num_obs - last_slope)) * sum(resids_adj.^2, 1);
      end
      		
      y_raw_ttl   (num_obs_ttl + 1:num_obs_ttl + num_obs) = y_raw_mth   (1:num_obs);
      y_adj_ttl   (num_obs_ttl + 1:num_obs_ttl + num_obs) = y_adj_mth   (1:num_obs);
      y_BA_adj_ttl(num_obs_ttl + 1:num_obs_ttl + num_obs) = y_BA_adj_mth(1:num_obs);
      y_hat_ttl   (num_obs_ttl + 1:num_obs_ttl + num_obs) = y_hat_mth   (1:num_obs);

      num_obs_ttl = num_obs_ttl + num_obs;
    end


        glbl_var_e(regr) = SSE / num_obs_ttl;
    glbl_var_raw_y(regr) = var(y_raw_ttl(1:num_obs_ttl));

    if style_passed == c.not_demeaned
      glbl_var_BA_adj_y(regr) = glbl_var_raw_y(regr);  
    else
      glbl_var_BA_adj_y(regr) = sum(y_BA_adj_ttl(1:num_obs_ttl).^2) ./ num_obs_ttl;  % demeaned monthly, so can sum squared values, adjusted for overfitting bias above divide by num_obs_ttl
    end

    glbl_r2_ttl_raw(regr) = 1 - glbl_var_e(regr)        / glbl_var_raw_y(regr);
    glbl_r2_dmn_raw(regr) = 1 - glbl_var_BA_adj_y(regr) / glbl_var_raw_y(regr);     % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
    glbl_r2_slps_raw_y  (regr) = glbl_r2_ttl_raw(regr) - glbl_r2_dmn_raw(regr);          % fraction of raw variance explained by x = R2 total - R2 demean
    glbl_r2_slps_dmn_y  (regr) = 1 - glbl_var_e(regr)        / glbl_var_BA_adj_y(regr);  % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y
      
       ave_r2(regr) = mean(r2_mth(frst_form_mth(regr_base):last_form_mth),    'omitnan');
    ave_var_e(regr) = mean(var_e_mth(frst_form_mth(regr_base):last_form_mth), 'omitnan');

    if 1 < regr
      ave_cor_PR_dP(regr - 1) = mean(cor_PR_dP(frst_form_mth(regr):last_form_mth), 'omitnan');
    end

    ave_cov_PR_LHS(regr) = mean(cov_PR_LHS(frst_form_mth(regr):last_form_mth), 'omitnan');
    
    ave_resid(:, regr) = mean(resid_mth, 'omitnan')';
             std_resid =  std(resid_mth, 'omitnan')';

    for area = 1:c.num_areas_both
      t_resid(area, regr) = ave_resid(area, regr) * sqrt(resid_cnt(area)) / std_resid(area);
    end
    
    ave_SD_x_adj(1:num_slopes, regr) = mean(SD_x_adj_mth(frst_form_mth(regr_base):last_form_mth, 1:num_slopes), 'omitnan');
     SD_SD_x_adj(1:num_slopes, regr) =  std(SD_x_adj_mth(frst_form_mth(regr_base):last_form_mth, 1:num_slopes), 'omitnan');
     
    if style_passed == c.not_demeaned
      [coef_stats(:, 1:num_slopes + 1, regr), autos(:, 1:num_slopes + 1, regr)] = summarize_FM_coefs(coef_mth(frst_form_mth(regr_base):last_form_mth, 1:num_slopes + 1), LHS_per, c, d, p);   % (3 + num_ts + num_autos) x num_coefs = (ave coef, se_smpl, se_LHS_1, ts, autos) x num_coefs
 
    else     % c.demean_CS_wo_int
      [coef_stats(:, 2:num_slopes + 1, regr), autos(:, 2:num_slopes + 1, regr)] = summarize_FM_coefs(coef_mth(frst_form_mth(regr_base):last_form_mth, 1:num_slopes), LHS_per, c, d, p);   % (3 + num_ts + num_autos) x num_coefs = (ave coef, se_smpl, se_LHS_1, ts, autos) x num_coefs
    end
  end

  FM_tmp.coef_stats          = coef_stats;
  FM_tmp.autos               = autos;
  FM_tmp.glbl_var_e          = glbl_var_e;
  FM_tmp.glbl_var_raw_y      = glbl_var_raw_y;
  FM_tmp.glbl_var_BA_adj_y   = glbl_var_BA_adj_y;
  FM_tmp.glbl_r2_ttl_raw     = glbl_r2_ttl_raw;
  FM_tmp.glbl_r2_dmn_raw     = glbl_r2_dmn_raw;
  FM_tmp.glbl_r2_slps_raw_y  = glbl_r2_slps_raw_y;
  FM_tmp.glbl_r2_slps_dmn_y  = glbl_r2_slps_dmn_y;

  FM_tmp.glbl_cov_pnl_wght_coef = NaN;
  FM_tmp.glbl_cor_pnl_wght_coef = NaN;
  FM_tmp.glbl_cor_var_x_coef    = NaN;
  
  FM_tmp.ave_r2       = ave_r2;
  FM_tmp.ave_var_e    = ave_var_e;
  FM_tmp.ave_SD_x_adj = ave_SD_x_adj;
  FM_tmp.SD_SD_x_adj  = SD_SD_x_adj;
  FM_tmp.cor_PR_dP    = ave_cor_PR_dP;           %       ave_cor_PR_dP(regr - 1) = mean(cor_pr_dP(frst_form_mth(regr):last_form_mth));
  FM_tmp.cov_PR_LHS   = ave_cov_PR_LHS;          %       ave_cov_PR_dP(regr) = mean(cor_pr_dP(frst_form_mth(regr):last_form_mth));

  FM_tmp.FE_r2        = (SSE_FE / num_obs_ttl) / glbl_var_raw_y(1);
  FM_tmp.ave_resid    = ave_resid;
  FM_tmp.t_resid      = t_resid;
end

% OK
%
function [coef, r2, var_e, resids_adj, y_hat] = cmpt_FM_coefs_and_r2(x_adj, y_raw, y_adj, source, nc, style_passed, c, d, p)

  num_obs   = size(x_adj, 1);
  num_coefs = size(x_adj, 2);
  
  % Compute coefficients
  %
  coef = (x_adj' * x_adj) \ x_adj' * y_adj;
  %Kx1
  
  coef = coef';
  %1xK


  % The TS average for each area may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y = var(y_raw);


  % Calculate residuals and adjust for overfitting if demeaned
  %
  y_hat  = x_adj * coef';
  resids = y_adj - y_hat;

  resids_adj = resids;

  % Calculate residual variances and r2
  %
  if style_passed == c.demean_CS_wo_int
    var_e = sum(resids_adj.^2, 1) / (num_obs - num_coefs - 1);
  else
    var_e = sum(resids_adj.^2, 1) / (num_obs - num_coefs);
  end
  r2 = 1.0 - (var_e ./ var_y);          % 1 x 1

end

% *******************************************************************************************************
%                                                                                                       *
%                                                End FM                                                 *
%                                                                                                       *
% *******************************************************************************************************

% *******************************************************************************************************
%                                                                                                       *
%                                          Panel Regressions                                            *
%                                                                                                       *
% *******************************************************************************************************
function Panel = run_Panel_dR_on_PR_and_one_dP(style_passed, c, d, p)

  % Runs Panel regressions of dR on PR and one dP. Called by Main.
  % "One" means for each combination of gap and dR, we run just PR, then PR and (one by one) the 1, 3, 6, and 12 mth change in dP to explain dR.
  % Calls run_set_of_Panel_dR_on_PR_and_one_dP for each combination of gap, LHS_per, and RHS_per. 

  % This function runs Panel regressions of dR on three variables -- the intercept, P/R, and dP.
  % dP might be the 1, 3, 6, or 12 month change in price ending in the same month P/R is observed.
  %
  % style_passed
  %
  % c.not_demeaned                   Panel    SE2  = {sum c {sum t e2ct}} / (N - K)
  %								  
  % c.demean_CS_w_int or _wo_int     Panel    SE2  = {sum t [Nt/(Nt - 1)] {sum c e2ct}} / (N - K)      
  %								  
  % c.demean_TS_w_int or _wo_int     Panel    SE2  = {sum c [Nc/(Nc - 1)] {sum t e2ct}} / (N - K)
  %
  % c.demean_both_CS_1st or _TS_1st  Panel    SE2  = {sum c {sum t [Nc/(Nc - 1)][Nt/(Nt - 1)]e2ct}} / (N - K)
  %
  % Demean_both makes sense only for panel...can include overall intercept
  %        TS 1st means set average TS value for each area  to zero, then set average CS value for each month to zero (TS average no longer zero)
  %        CS 1st means set average CS value for each month to zero, then set average TS value for each area  to zero (CS average no longer zero)
  %
  % c.not_demeaned          = 1; x
  % c.demean_CS_wo_int      = 3; x       % This makes sense only for TS and panel...without intercept
  %

  if ~x_in_y(style_passed, [c.not_demeaned, c.demean_CS_wo_int])
    disp('style_passed to run_Panel_dR_on_PR_and_one_dP does not make sense')
    disp(c.demean_label(style_passed))
    pause
  end
  
  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.Panel, style_passed, c, d, p);
							 
  num_coef_stats = get_num_Panel_coef_stats;  % 11

  coef_stats_label  = load_Panel_coef_stats_label;
  
  Panel.coef_stats         = zeros(num_coef_stats, max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);   % need max_num_coefs because run_Panel returns results for intercept
  Panel.ave_SD_x_adj       =                zeros(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);   
  Panel.SD_SD_x_adj        =                zeros(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);   

  Panel.glbl_r2_ttl_raw    =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % 1 - variance of bias adjusted residuals relative to variance of raw y
  Panel.glbl_r2_dmn_raw    =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y
  Panel.glbl_r2_slps_raw_y =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by x = R2 total - R2 demean
  Panel.glbl_r2_slps_dmn_y =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y

  Panel.glbl_var_e         =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  Panel.glbl_var_raw_y     =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  Panel.glbl_var_BA_adj_y  =                                zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  % Cycle over gaps
  %
  for LHS_RHS_var_indx = 1:num_LHS_RHS_combos
    LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
    RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);
      
    for LHS_per_indx = 1:num_LHS_pers
      LHS_per = p.LHS_per(LHS_per_indx);

      Panel_tmp = run_set_of_Panel_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p);

               Panel.coef_stats(:, :, :, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.coef_stats;
                Panel.glbl_r2_ttl_raw(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_r2_ttl_raw;
                Panel.glbl_r2_dmn_raw(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_r2_dmn_raw;
                  Panel.glbl_r2_slps_raw_y(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_r2_slps_raw_y;
                  Panel.glbl_r2_slps_dmn_y(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_r2_slps_dmn_y;
 
                     Panel.glbl_var_e(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_var_e;
                 Panel.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_var_raw_y;
              Panel.glbl_var_BA_adj_y(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.glbl_var_BA_adj_y;

      Panel.glbl_cov_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = NaN;
      Panel.glbl_cor_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = NaN;
         Panel.glbl_cor_var_x_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = NaN;
        
                Panel.ave_SD_x_adj(:, :, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.ave_SD_x_adj;
                 Panel.SD_SD_x_adj(:, :, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.SD_SD_x_adj;
                      Panel.cor_PR_dP(:, LHS_per_indx, LHS_RHS_var_indx) = Panel_tmp.cor_PR_dP;
    end
  end

  %  If x and y demeaned...no intercept in panel...intercept inferred from means and slopes
  %     coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both
  %  else
  %     coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both
  %
  % drop column for intercept
  %
  Panel.coef_stats = Panel.coef_stats(:, 2:end, :, :, :, :, :, :);
end


% OK
%
function Panel = run_set_of_Panel_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p)

  % Set here means all regressions for combination of gap and LHS_per, we run just PR, then PR and (one by one) the 1, 3, 6, and 12 mth change in dP to explain dR for all months.
  % Computes coefs and r2s directly in the function.

  if ~x_in_y(style_passed, [c.not_demeaned, c.demean_CS_wo_int])
    disp('style_passed to run_set_of_Panel_dR_on_PR_and_one_dP does not make sense')
    disp(style_passed)
    pause
  end

  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.FM, style_passed, c, d, p);


  % coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both
  %
  num_coef_stats = get_num_Panel_coef_stats;  %  11
  
  resids       =              NaN(c.max_TS_CS_obs,  num_regrs);
  nt           =             zeros(c.num_mths_ttl,  num_regrs);
  nc           =           zeros(c.num_areas_both,  num_regrs);
  y_raw        =              NaN(c.max_TS_CS_obs,  num_regrs);
  source       =              NaN(c.max_TS_CS_obs,  num_regrs);

  autos        = NaN(num_autos,      max_num_slopes, num_regrs);

  % Global results
  %
  Panel.coef_stats        = NaN(num_coef_stats, max_num_coefs,  num_regrs);
  Panel.glbl_r2_ttl_raw   =                                 NaN(num_regrs, 1);
  Panel.glbl_r2_dmn_raw   =                                 NaN(num_regrs, 1);
  Panel.glbl_r2_slps_raw_y     =                                 NaN(num_regrs, 1);
  Panel.glbl_r2_slps_dmn_y     =                                 NaN(num_regrs, 1);

  Panel.glbl_var_e        =                                 NaN(num_regrs, 1);
  Panel.glbl_var_raw_y    =                                 NaN(num_regrs, 1);
  Panel.glbl_var_BA_adj_y =                                 NaN(num_regrs, 1);
										
  Panel.ave_SD_x_adj      =                 NaN(max_num_slopes, num_regrs);
  Panel.SD_SD_x_adj       =                 NaN(max_num_slopes, num_regrs);
					  
  Panel.cor_PR_dP         = NaN(num_RHS_pers, 1);
  
  
  % Working arrays
  %

  y_raw    = NaN(c.max_TS_CS_obs, 1);
  y_adj    = NaN(c.max_TS_CS_obs, 1);
  y_BA_adj = NaN(c.max_TS_CS_obs, 1);

  x_raw    = NaN(c.max_TS_CS_obs, max_num_coefs);
  x_adj    = NaN(c.max_TS_CS_obs, max_num_coefs);

  source   = NaN(c.max_TS_CS_obs, 1);

  if style_passed == c.not_demeaned
    x_adj(:, 1) = 1;
  end


  % determine and possibly write frst and last form mth
  %
  % frst_form_mth = NaN(num_regrs, 1);  
  % last_form_mth = NaN;                       
  % frst_form_mth_area = NaN(c.num_areas_both, num_regrs);    % num_regrs...1 is for simple P/R
  % last_form_mth_area = NaN(c.num_areas_both, 1);            

  [frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, nt, nc] = determine_and_wrt_frst_and_last_form_mth(LHS_per, LHS_var, RHS_var, c.dont_wrt_frst_last_mth, c, d, p);
  

  % Select LHS_chg data based on LHS_per
  %
  LHS_raw = load_mthly_qtrly_semi_or_annual_chg(LHS_var, LHS_per, 0, c, d, p);  % The two omitted results are for subtracting the overall mean for the month (in the TS regression)
  

  % Load vector of PR or RP data
  %
  PR_raw = load_full_PR_vector_regr(0, c, d, p);


  % Run Panel regression with only PR if regr = 1 or PR and dP if regr > 1
  %
  for regr = 1:num_regrs
    
    % Setup RHS variable
    %
    if regr == 1
      last_slope = frst_slope;
      num_slopes = 1;
      regr_base  = regr;
      
    else
      last_slope = frst_slope + 1;
      num_slopes = 2;
      regr_base  = regr;
      
    end
    
    if regr ~= 1
      RHS_per = p.RHS_per(regr_base - 1);
      RHS_raw_lag  = load_mthly_qtrly_semi_or_annual_chg(RHS_var, RHS_per, 0, c, d, p);  % This loads RHS_var for all areas
      RHS_raw_lead = load_mthly_qtrly_semi_or_annual_chg(RHS_var, LHS_per, 0, c, d, p);  % This loads RHS_var for all areas
      
      if style_passed == c.not_demeaned
        RHS_adj_lag = RHS_raw_lag;
        
      else
        RHS_adj_lag  = Panel_demean_x(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, RHS_raw_lag, style_passed, c, d, p);
        RHS_adj_lead = Panel_demean_y(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, RHS_raw_lead, style_passed, c, d, p);
      end
    end

      
    % Load demeaned values if needed...Note -- PR_adj(mth) explains LHS_adj(mth + LHS_per)
    %
    PR_adj  = Panel_demean_x(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, PR_raw, style_passed, c, d, p);
    LHS_adj = Panel_demean_y(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, LHS_raw, style_passed, c, d, p);
  

    num_obs  = 0;
    frst_mth = frst_form_mth(regr_base);
    last_mth = last_form_mth;    
    for mth = frst_form_mth(regr_base):last_form_mth
      mth_ttl = 0;
      for area = 1:c.num_areas_both
        if (0 < min(d.frst_ob(area, :))) && (frst_form_mth_area(area, regr_base) <= mth) && (mth <= last_form_mth_area(area))
          num_obs                      = num_obs + 1;
          mth_ttl                      = mth_ttl + 1;
          y_raw(num_obs)               = LHS_raw(mth + LHS_per, area);
          y_adj(num_obs)               = LHS_adj(mth + LHS_per, area);

          if style_passed == c.not_demeaned
            y_BA_adj(num_obs)          = LHS_adj(mth + LHS_per, area);
          elseif style_passed == c.demean_CS_wo_int
            y_BA_adj(num_obs)          = LHS_adj(mth + LHS_per, area) * sqrt(nt(mth, regr_base) / (nt(mth, regr_base) - 1));
          end
          
          x_raw(num_obs, frst_slope) = PR_raw(mth, area);
          x_adj(num_obs, frst_slope) = PR_adj(mth, area);
          
          if regr ~= 1
            x_raw(num_obs, frst_slope + 1) = RHS_raw_lag(mth, area);
            x_adj(num_obs, frst_slope + 1) = RHS_adj_lag(mth, area);
          end

          source(num_obs) = area;
        end
      end
    end

    if 1 < regr 
      Panel.cor_PR_dP(regr - 1) = corr(x_adj(1:num_obs, frst_slope), x_adj(1:num_obs, frst_slope + 1));
    end
     
    Panel_shrt = cmpt_Panel_regr(x_raw(1:num_obs, 1:last_slope), x_adj(1:num_obs, 1:last_slope), y_raw(1:num_obs), y_adj(1:num_obs), y_BA_adj(1:num_obs), source(1:num_obs), LHS_per, nt(frst_mth:last_mth, regr_base), nc(:, regr_base), style_passed, c, d, p);

    Panel.coef_stats(:, 1:num_slopes + 1, regr) = Panel_shrt.coef_stats(1:num_coef_stats, 1:num_slopes + 1);   % num_slopes + 1 because intercept

    Panel.glbl_var_e       (regr) = Panel_shrt.glbl_var_e;       
    Panel.glbl_var_raw_y   (regr) = Panel_shrt.glbl_var_raw_y;   
    Panel.glbl_var_BA_adj_y(regr) = Panel_shrt.glbl_var_BA_adj_y;
						   		
    Panel.glbl_r2_ttl_raw  (regr) = Panel_shrt.glbl_r2_ttl_raw;  
    Panel.glbl_r2_dmn_raw  (regr) = Panel_shrt.glbl_r2_dmn_raw;    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
    Panel.glbl_r2_slps_raw_y    (regr) = Panel_shrt.glbl_r2_slps_raw_y;      % fraction of raw variance explained by x = R2 total - R2 demean
    Panel.glbl_r2_slps_dmn_y    (regr) = Panel_shrt.glbl_r2_slps_dmn_y;      % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y
								
    Panel.ave_SD_x_adj(1:num_slopes, regr) = Panel_shrt.ave_SD_x_adj;     
    Panel.SD_SD_x_adj (1:num_slopes, regr) = Panel_shrt.SD_SD_x_adj;      

  end
end


function Panel = cmpt_Panel_regr(x_raw, x_adj, y_raw, y_adj, y_BA_adj, source, LHS_per, nt, nc, style_passed, c, d, p)

       % Kx1  KxK        KxK            KxK       KxK         1    1

  %     coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both
  %
  %     col 1 summarizes monthly intercepts
  %         coef_stats(1,  1) = coef_stats_shrt(1, 1);  % ave int
  %         coef_stats(4,  1) = coef_stats_shrt(3, 1);  % SE ave int
  %         coef_stats(11, 1) = coef_stats_shrt(4, 1);  % t-simple
  %         coef_stats(12, 1) = coef_stats_shrt(5, 1);  % t LHS_per-1
  %
  %     coef_stats = [P'; SE_P; t_stats];                      % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time, area, both;  t-stats: iid, hetero, clstrd time, area, both
  %
  %     If intercept is included in the panel the estimated intercept is summarized in col
  %         coef_stats(1,  1) = coef_stats_shrt(1, 1);  % coefficient
  %         coef_stats(4,  1) = coef_stats_shrt(3, 1);  % SE ave int
  %         coef_stats(11, 1) = coef_stats_shrt(4, 1);  % t-simple
  %         coef_stats(12, 1) = coef_stats_shrt(5, 1);  % t LHS_per-1
  %
  % For convenience, this discussion assumes data are clustered by period, with period dummies.
  % It does not matter, however, whether grouping is by periods, firms, areas, etc.
  %
  % The White standard errors control for heteroscedasticity.
  % The time clustered standard errors assume only the contemporaneous residuals are correlated.
  % The area clustered standard errors assume each area's residuals are correlated through time, but there is no cross-area correlation.
  % The both clustered standard errors assume each area's residuals are correlated through time and there is contemporaneous cross-area correlation.
  % The panel can be unbalanced, with different number of observations across periods.
  % x and y must be grouped by period
  % There are T periods and n(t) is the number of observations in period t.
  % The first n(1) observations are for period 1, the next n(2) observations are for period 2, etc.
  % sum(n) must be size(x, 1) = size(y, 1).
  %
  % ave_SD_x_adj and SD_SD_x_adj are the average and standard deviation of the monthly CS standard deviations of x_adj (demeaned x)
  %

  num_mths  = size(nt,    1);
  num_obs   = size(x_adj, 1);
  num_coefs = size(x_adj, 2);

  resids    = NaN(c.max_TS_CS_obs, 1);
  intcpts   = NaN(num_mths, 1);  

  if sum(nt) ~= num_obs
    disp('Sum of nt must equal rows in x and y')
    'sum(nt) ', sum(nt)
    num_obs
    LHS_per
    pause
  end
  
  [ave_SD_x_adj, SD_SD_x_adj] = summarize_CS_SD_of_x_adj(x_adj, nt, style_passed, c, d, p);  % Column 1 for ave_SD and SD_SD is for frst_slope  num_slopes x 1


  % compute coefs and residuals
  %
  xx_adj = x_adj' * x_adj;                                      % num_coefs x num_coefs
  P      = xx_adj \ x_adj' * y_adj;                             % num_coefs x 1

  resids = y_adj - x_adj * P;
  
  var_raw_y = var(y_raw);

  if style_passed == c.not_demeaned
    var_BA_adj_y = var_raw_y;
  else
    var_BA_adj_y = sum(y_BA_adj.^2) ./ num_obs;  % demeaned monthly, so can sum squared values, adjusted for overfitting bias above divide by num_obs_ttl
  end

  
  % compute var_e, r2, and residuals adjusted for demeaning

  [var_e, r2, resids_adj] = cmpt_Panel_var_e_r2_resids_adj(x_adj, resids, nt, nc, source, var_raw_y, style_passed, c, d, p);
  
  
  % Compute SEs of Panel coefs
  %                                                                                                       2                         3                        4                                           5                              6
  SE_P = cmpt_Panel_SEs_of_coefs(x_adj, var_e, resids_adj, nt, LHS_per, source, c, d, p);   % SE_P = sqrt[(diag(iid_cov_P)'); (diag(cov_P_white)'); (diag(clstrd_cov_P_time) or generalized HH'); (diag(clstrd_cov_P_area)'); (diag(clstrd_cov_P_both)')];   % 5 x K
  
  t_stats = P' ./ SE_P;                                                 % 5 x K 
                                                                        %                    1           2     3                  4                  5     6               7     8                   9                  10   11
  coef_stats = [P'; SE_P; t_stats];                                     % 11 x K    Rows... Coef;  SEs: iid, hetero, clstrd time or generalized HH, area, both;  t-stats: iid, hetero, clstrd time or generalized HH, area, both


  % add summary stats for implied intercepts if not included in x
  %
  if x_in_y(style_passed, [c.demean_CS_wo_int, c.demean_TS_wo_int, c.demean_both_TS_1st_wo, c.demean_both_CS_1st_wo, c.demean_TS_wo_int, c.demean_both_CS_1st_wo, c.demean_both_full])

    coef_stats = add_sum_stats_for_implied_intcept_Panel(y_raw, x_raw, nt, source, LHS_per, P, style_passed, coef_stats, c, d, p);

  end
  
  Panel.coef_stats        = coef_stats;
  Panel.glbl_var_e        = var_e;
  Panel.glbl_var_raw_y    = var_raw_y;
  Panel.glbl_var_BA_adj_y = var_BA_adj_y;

  Panel.glbl_r2_ttl_raw   = 1 - var_e        / var_raw_y;
  Panel.glbl_r2_dmn_raw   = 1 - var_BA_adj_y / var_raw_y;                   % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
  Panel.glbl_r2_slps_raw_y     = Panel.glbl_r2_ttl_raw - Panel.glbl_r2_dmn_raw;  % fraction of raw variance explained by x = R2 total - R2 demean
  Panel.glbl_r2_slps_dmn_y     = 1 - var_e        / var_BA_adj_y;                % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y

  Panel.ave_SD_x_adj      = ave_SD_x_adj;
  Panel.SD_SD_x_adj       = SD_SD_x_adj;
end

function  [var_e, r2, resids_adj] = cmpt_Panel_var_e_r2_resids_adj(x_adj, resids, nt, nc, source, var_raw_y, style_passed, c, d, p)

  num_mths  = size(nt,    1);
  num_obs   = size(x_adj, 1);
  num_coefs = size(x_adj, 2);

  resids_adj = NaN(size(resids));
  
  % Adjust residuals for demeaning and compute var_e and R2
  %
  if style_passed == c.not_demeaned
    resids_adj = resids;
    sse        = sum(resids.^2, 1);
    
  elseif style_passed == c.demean_CS_wo_int
    sse        = 0;
    obs_end    = 0;
    for t = 1:num_mths
      obs_beg = obs_end + 1;           
      obs_end = obs_end + nt(t);
      sse     = sse + (nt(t)/(nt(t) - 1)) * sum(resids(obs_beg:obs_end).^2, 1);

      resids_adj(obs_beg:obs_end) = sqrt((nt(t)/(nt(t) - 1))) * resids(obs_beg:obs_end);
    end
  end

  var_e = sse / (num_obs - num_coefs);
  
  r2    = 1 - (var_e ./ var_raw_y);

end


function SE_P = cmpt_Panel_SEs_of_coefs(x_adj, var_e, resids_adj, nt, LHS_per, source, c, d, p)

  % SE_P = sqrt[(diag(iid_cov_P)'); (diag(cov_P_white)'); (diag(clstrd_cov_P_time)'); (diag(clstrd_cov_P_area)'); (diag(clstrd_cov_P_both)')];   % 5 x K
  
  num_mths  = size(nt,    1);
  num_obs   = size(x_adj, 1);
  num_coefs = size(x_adj, 2);
  
  n_beg = zeros(num_mths, 1);
  n_beg(1) = 1;
  for t = 2:num_mths
    n_beg(t) = n_beg(t-1) + nt(t-1);
  end
  n_end = [n_beg(2:t) - 1; num_obs];



  xx_adj    = x_adj' * x_adj;                                     % num_coefs x num_coefs

  % Compute SEs of coefs assuming iid resids
  %
  inv_xx_adj = inv(xx_adj);
  iid_cov_P  = var_e * inv_xx_adj;
  
  
  % Compute White heteroscedasticity-consistent SEs
  %
  x_omega_x_white = x_adj(1:num_obs, :)' * (x_adj(1:num_obs, :) .* (resids_adj(1:num_obs, 1).^2));
  cov_P_white = inv_xx_adj * x_omega_x_white * inv_xx_adj;    % num_coefs x num_coefs


  % Compute time-clustered SEs = Inv(x_adj' * x_adj) * x_omega_x * Inv(x_adj' * x_adj)
  %
  x_omega_x_time = zeros(num_coefs, num_coefs);
  if p.panel_time_clstrd_just_contemp
    obs_end   = 0;
    for t = 1:num_mths
      obs_beg = obs_end + 1;
      obs_end = obs_end + nt(t);

      x_omega_x_time = x_omega_x_time + x_adj(obs_beg:obs_end, :)' * resids_adj(obs_beg:obs_end, 1) * resids_adj(obs_beg:obs_end, 1)' * x_adj(obs_beg:obs_end, :);   % num_coefs x num_coefs
    end


  % Generalized HH
  %
  else
    for t1 = 1:num_mths
      for t2 = max(1, t1 - LHS_per + 1):min(num_mths, t1 + LHS_per - 1)
        x_omega_x_time = x_omega_x_time + x_adj(n_beg(t1):n_end(t1), :)' * resids_adj(n_beg(t1):n_end(t1), :) * resids_adj(n_beg(t2):n_end(t2), :)' * x_adj(n_beg(t2):n_end(t2), :);
      end
    end      
  end
  
  clstrd_cov_P_time = inv_xx_adj * x_omega_x_time * inv_xx_adj;    % num_coefs x num_coefs...this is genrlzd_HH_cov_TS in TS results
      


  % Compute area-clustered SEs = Inv(x_adj' * x_adj) * x_omega_x * Inv(x_adj' * x_adj)
  %
  x_omega_x_area = zeros(num_coefs, num_coefs);
  area_obs       = zeros(c.num_areas_both, 1);
  area_x_adj     = NaN(num_obs, num_coefs);
  area_resid_adj = NaN(num_obs, num_coefs);
  
  for area = 1:c.num_areas_both
    for obs = 1:num_obs
      if source(obs, 1) == area
        area_obs(area) = area_obs(area) + 1;
        
            area_x_adj(area_obs(area), :) = x_adj(obs, :);
        area_resid_adj(area_obs(area), 1) = resids_adj(obs, 1);
      end
    end
      
    x_omega_x_area = x_omega_x_area + area_x_adj(1:area_obs(area), :)' * area_resid_adj(1:area_obs(area), 1) * area_resid_adj(1:area_obs(area), 1)' * area_x_adj(1:area_obs(area), :);   % num_coefs x num_coefs
  end

  clstrd_cov_P_area = inv_xx_adj * x_omega_x_area * inv_xx_adj;    % num_coefs x num_coefs
  
  x_omega_x_both = x_omega_x_time + x_omega_x_area - x_omega_x_white;
  
  clstrd_cov_P_both = inv_xx_adj * x_omega_x_both * inv_xx_adj;    % num_coefs x num_coefs

  SE_P = sqrt([diag(iid_cov_P)'; diag(cov_P_white)'; diag(clstrd_cov_P_time)'; diag(clstrd_cov_P_area)'; diag(clstrd_cov_P_both)']);   % 5 x K

end

%
% *****************************************************************************************
%                                                                                         *
%                                  End Panel Regressions                                  *
%                                                                                         *
% *****************************************************************************************
%
% ****************************************************************************************
%                                                                                        *
%                                    TS Regressions                                      *
%                                                                                        *
% ****************************************************************************************
%
function TS = run_TS_dR_on_PR_and_one_dP(style_passed, c, d, p)

  % Demeaning here means subtract the CS average of x or y...Since this is a TS regression, that does not mean the averages of x and y are zero. 
  % 
  % style_passed
  %
  % c.not_demeaned                   TS       SE2c = {sum t e2ct} / (Nc - K)                        SE2 = sum c {Nc SE2c / N}
  %								  
  % c.demean_CS_w_int or _wo_int     TS       SE2c = {sum t [Nt/(Nt - 1)]e2ct} / (Nc - K)           SE2 = sum c {Nc SE2c / N}
  %								  

  % 
  % coef_stats    = NaN(p.num_coef_stats_simple, max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);     % FM, with areas like months...ave_coef, sd_coef, coef_SE, t-simple, last (t-1) empty

  % p.RHS_per = [1, 3, 6, 12];
  % p.LHS_per = [1, 3, 6, 12];

  if ~x_in_y(style_passed, [c.not_demeaned, c.demean_CS_w_int, c.demean_CS_wo_int, c.demean_both_full])
    disp('style_passed must be c.not_demeaned, c.demean_CS_w_int, c.demean_CS_wo_int, or c.demean_TS_wo_int in run and wrt_TS_dR_on_PR_and_one_dP')
    disp(c.demean_label(style_passed))
    pause
  end

  num_areas = c.num_areas_both;

  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.TS, style_passed, c, d, p);

  num_autos      = p.num_autos_resids;

  [num_coef_stats, num_area_coef_stats, num_ave_ts, num_area_ts] = get_num_TS_coef_stats(p);   % num_coef_stats is for the CS ave of multiple TS regressions, num_area_coef_stats is for one TS regression.

  PR_label       = load_PR_label(c, p);
  PR_label_shrt  = load_PR_label_shrt(c, p);
  
  if p.frst_last_form_mth_same_or_diff == c.same
    same_or_diff_label = 'Same';
  else
    same_or_diff_label = 'Diff';
  end


  % Global (cross-area) results
  %
  TS.coef_stats        = NaN(num_coef_stats, max_num_coefs,  num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.glbl_r2_ttl_raw   =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % 1 - variance of bias adjusted residuals relative to variance of raw y
  TS.glbl_r2_dmn_raw   =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y
  TS.glbl_r2_slps_raw_y =                              zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of raw variance explained by x = R2 total - R2 demean
  TS.glbl_r2_slps_dmn_y =                              zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);    % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y
  TS.glbl_var_e        =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.glbl_var_raw_y    =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.glbl_var_BA_adj_y =                               zeros(num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  TS.ave_r2_raw_y      =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.ave_r2_dmn_y      =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.ave_var_e         =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.ave_SD_x_adj      =                 NaN(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.SD_SD_x_adj       =                 NaN(max_num_slopes, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.ave_resid_cor     =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.cor_PR_dP         =                              NaN(num_RHS_pers, num_LHS_pers);
  TS.FE_r2             =                                            NaN(num_LHS_pers, num_LHS_RHS_combos);

  TS.glbl_cov_pnl_wght_coef =             NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.glbl_cor_pnl_wght_coef =             NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  TS.glbl_cor_var_x_coef    =             NaN(max_num_coefs, num_regrs, num_LHS_pers, num_LHS_RHS_combos);



  % Area results
  %
  TS.coef_stats_area = NaN(num_area_coef_stats, max_num_coefs,  num_regrs, c.num_areas_both, num_LHS_pers, num_LHS_RHS_combos);
  TS.autos_area      =                           NaN(num_autos, num_regrs, c.num_areas_both, num_LHS_pers, num_LHS_RHS_combos);
  TS.r2_area         =                                      NaN(num_regrs, c.num_areas_both, num_LHS_pers, num_LHS_RHS_combos);
  TS.var_e_area      =                                      NaN(num_regrs, c.num_areas_both, num_LHS_pers, num_LHS_RHS_combos);
  TS.SD_x_adj_area   =                      NaN(max_num_slopes, num_regrs, c.num_areas_both, num_LHS_pers, num_LHS_RHS_combos);
  TS.cor_PR_dP_area  =                                   NaN(num_RHS_pers, c.num_areas_both, num_LHS_pers);
 
  % Cycle over LHS_RHS_combos and LHS periods
  %
  for LHS_RHS_var_indx = 1:num_LHS_RHS_combos
    LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
    RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);

    for LHS_per_indx = 1:num_LHS_pers
      LHS_per = p.LHS_per(LHS_per_indx);

      TS_tmp = run_set_of_TS_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p);

      % Global results
      %
      TS.coef_stats(:, :,  :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.coef_stats;

      TS.glbl_r2_ttl_raw   (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_r2_ttl_raw;
      TS.glbl_r2_dmn_raw   (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_r2_dmn_raw;
      TS.glbl_r2_slps_raw_y(:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_r2_slps_raw_y;
      TS.glbl_r2_slps_dmn_y(:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_r2_slps_dmn_y;

      TS.glbl_var_e       (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_var_e;
      TS.glbl_var_raw_y   (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_var_raw_y;
      TS.glbl_var_BA_adj_y(:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_var_BA_adj_y;

      TS.ave_SD_x_adj (:,  :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.ave_SD_x_adj;
      TS.SD_SD_x_adj  (:,  :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.SD_SD_x_adj;

      TS.ave_r2_raw_y     (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.ave_r2_raw_y;
      TS.ave_r2_dmn_y     (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.ave_r2_dmn_y;
      TS.ave_var_e        (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.ave_var_e;
      TS.ave_resid_cor    (:, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.ave_resid_cor;
      TS.cor_PR_dP        (:, LHS_per_indx)                   = TS_tmp.cor_PR_dP;
      
      TS.glbl_cov_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_cov_pnl_wght_coef;
      TS.glbl_cor_pnl_wght_coef(:, :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_cor_pnl_wght_coef;
      TS.glbl_cor_var_x_coef   (:, :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.glbl_cor_var_x_coef;
          
      % Area results
      %
      TS.coef_stats_area(:, :, :, :, LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.coef_stats_area;
      TS.autos_area     (:, :, :,    LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.autos_area;
      TS.r2_raw_y_area  (:, :,       LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.r2_raw_y_area;
      TS.r2_dmn_y_area  (:, :,       LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.r2_dmn_y_area;
      TS.var_e_area     (:, :,       LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.var_e_area;
      TS.SD_x_adj_area  (:, :, :,    LHS_per_indx, LHS_RHS_var_indx) = TS_tmp.SD_x_adj_area;
      TS.cor_PR_dP_area    (:, :,    LHS_per_indx)                   = TS_tmp.cor_PR_dP_area;

    end
  end
end

% OK
%
function wrt_ave_TS_regr_results(LHS_RHS_var_indx, TS, num_coef_stats, coef_stats_label, style_passed, c, d, p)   

  FM_Panel_or_TS_label = 'TS';
  
  LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
  RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);

  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.TS, style_passed, c, d, p);
							
  if LHS_var == RHS_var
    num_regrs = num_RHS_pers + 1;
  end
							
  RHS_per_list   = [0, p.RHS_per];
  LHS_per_list   = p.LHS_per;
  max_RHS_per    = max(RHS_per_list);
  
  PR_label_shrt  = load_PR_label_shrt(c, p);

  LHS_dZ_label   = load_dP_or_dR_label(LHS_var, c);
  RHS_dZ_label   = load_dP_or_dR_label(RHS_var, c);

  frst_form_mth  = NaN(1 + num_RHS_pers, 1);
  
  last_slope     = frst_slope + max_num_slopes - 1;
 
  out_format = [p.slope_format, p.se_format, p.se_format, p.t_stat_format, p.t_stat_format];
 
  % Cycle over LHS_chg periods
  %
  for LHS_per_indx = 1:num_LHS_pers
    LHS_per = LHS_per_list(LHS_per_indx);
                            
    if p.frst_last_form_mth_same_or_diff == c.same
      [frst_form_mth(1), last_form_mth] = cmpt_frst_and_last_formation_mths(0, LHS_var, RHS_var, LHS_per, max_RHS_per, c, d, p, c.PR_reqd);  % frst_form_mth scalar, last_form_mth scalar determined by LHS_per
      frst_form_mth = frst_form_mth(1);
    else
      [frst_form_mth(:), last_form_mth] = cmpt_frst_and_last_formation_mths(0, LHS_var, RHS_var, LHS_per, RHS_per_list, c, d, p, c.PR_reqd);  % frst_form_mth 1 x num_RHS_pers, last_form_mth scalar determined by LHS_per
    end
  
    fprintf(p.fid, '\n\nLHS Period = %d     %s on %s and %s   %s  (wrt_ave_TS_regr_results)', LHS_per, LHS_dZ_label, PR_label_shrt, RHS_dZ_label, c.demean_label(style_passed));
    
    if p.frst_last_form_mth_same_or_diff == c.same
      fprintf(p.fid, '\n                                    Same Start Months    %6d-%6d', d.idate(frst_form_mth(1)), d.idate(last_form_mth));
      fprintf(p.fid, '\n                    '); for regr = 1:num_RHS_pers; fprintf(p.fid, '      RHS Period =%2d', p.RHS_per(regr)); end
    else
      fprintf(p.fid, '\n\n                                          Different Start Months');
      fprintf(p.fid, '\n            %6d  ', d.idate(frst_form_mth(1))); for regr = 1:num_RHS_pers; fprintf(p.fid, '      RHS Period =%2d', p.RHS_per(regr)); end; for regr = 1:num_RHS_pers; fprintf(p.fid, '           RHS Period =%2d   ', p.RHS_per(regr)); end
      fprintf(p.fid, '\n             -%6d', d.idate(last_form_mth(1))); 
      for regr = 2:num_RHS_pers + 1;     fprintf(p.fid, '       %6d-%6d', d.idate(frst_form_mth(regr)), d.idate(last_form_mth)); end
      for regr = 2:num_RHS_pers + 1;     fprintf(p.fid, '           %6d-%6d    ', d.idate(frst_form_mth(regr)), d.idate(last_form_mth)); end
    end
              
    frst_coef = 1;
    fprintf(p.fid, '\n                %s ', PR_label_shrt); for regr = 1:num_RHS_pers; fprintf(p.fid, '         %s   %s-K', PR_label_shrt, RHS_dZ_label); end; for regr = 1:num_RHS_pers; fprintf(p.fid, '         %s   %s-K   %s+T', PR_label_shrt, RHS_dZ_label, RHS_dZ_label); end
    
    for coef_stats_indx = 1:num_coef_stats

      if out_format(coef_stats_indx) == c.f82
        fprintf(p.fid, '\n%s %6.2f     ', coef_stats_label(coef_stats_indx, 1), TS.coef_stats(coef_stats_indx, frst_slope, 1, LHS_per_indx, LHS_RHS_var_indx)); 
        for regr = 2:num_RHS_pers + 1;         fprintf(p.fid, '%8.2f',  TS.coef_stats(coef_stats_indx, frst_slope:frst_slope + 1, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end
        for regr = num_RHS_pers + 2:num_regrs; fprintf(p.fid, '%8.2f',  TS.coef_stats(coef_stats_indx, frst_slope:frst_slope + 2, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end


      elseif out_format(coef_stats_indx) == c.f83
        fprintf(p.fid, '\n%s %6.3f     ', coef_stats_label(coef_stats_indx, 1), TS.coef_stats(coef_stats_indx, frst_slope, 1, LHS_per_indx, LHS_RHS_var_indx)); 
        for regr = 2:num_RHS_pers + 1;         fprintf(p.fid, '%8.3f',  TS.coef_stats(coef_stats_indx, frst_slope:frst_slope + 1, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end
        for regr = num_RHS_pers + 2:num_regrs; fprintf(p.fid, '%8.3f',  TS.coef_stats(coef_stats_indx, frst_slope:frst_slope + 2, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end

      else
        disp('p.out_format must be c.f82 or c.f83 in wrt_ave_TS_regr_results')
        pause
      end
    end
		
    r2_shrt = squeeze(TS.glbl_r2_ttl_raw(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nR2 Ttl Raw   %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    r2_shrt = squeeze(TS.glbl_r2_dmn_raw(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nR2 Dmn Raw   %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    r2_shrt = squeeze(TS.glbl_r2_slps_raw_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nR2 Slps Raw  %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    r2_shrt = squeeze(TS.glbl_r2_slps_dmn_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nR2 Slps Dmn  %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    r2_shrt = squeeze(TS.ave_r2_raw_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nAve R2 Raw Y %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    r2_shrt = squeeze(TS.ave_r2_dmn_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nAve R2 Dmn Y %6.2f     ', r2_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', r2_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    var_e_shrt = squeeze(TS.glbl_var_e(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nGlobal Var e %6.2f     ', var_e_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', var_e_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', var_e_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    var_y_shrt = squeeze(TS.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nGlobal Var Y %6.2f     ', var_y_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', var_y_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', var_y_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    var_y_shrt = squeeze(TS.glbl_var_BA_adj_y(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nVar BA Adj Y %6.2f     ', var_y_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', var_y_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', var_y_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    SD_e_shrt = sqrt(squeeze(TS.glbl_var_e(:, LHS_per_indx, LHS_RHS_var_indx)));
    fprintf(p.fid, '\nGlobal SD e  %6.2f     ', SD_e_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', SD_e_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', SD_e_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    SD_y_shrt = sqrt(squeeze(TS.glbl_var_raw_y(:, LHS_per_indx, LHS_RHS_var_indx)));
    fprintf(p.fid, '\nGlobal SD Y  %6.2f     ', SD_y_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', SD_y_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', SD_y_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    SD_y_shrt = sqrt(squeeze(TS.glbl_var_BA_adj_y(:, LHS_per_indx, LHS_RHS_var_indx)));
    fprintf(p.fid, '\nGlobal SD Y  %6.2f     ', SD_y_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', SD_y_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', SD_y_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    var_e_shrt = squeeze(TS.ave_var_e(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nAve Var e    %6.2f     ', var_e_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f', var_e_shrt(regr)); fprintf(p.fid, '            '); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f', var_e_shrt(regr)); for coef = frst_slope + 1:last_slope; fprintf(p.fid, '        '); end; fprintf(p.fid, '    '); end

    fprintf(p.fid, '\nAve SD adj x %6.2f     ',  TS.ave_SD_x_adj(1, LHS_per_indx, LHS_RHS_var_indx)); 
    for regr = 2:num_regrs_base; fprintf(p.fid, '%8.2f', TS.ave_SD_x_adj(1:2, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end
    for regr = num_regrs_base + 1:num_regrs; fprintf(p.fid, '%8.2f', TS.ave_SD_x_adj(1:3, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end

    fprintf(p.fid, '\nSD  SD adj X %6.2f     ',  TS.SD_SD_x_adj(1, LHS_per_indx, LHS_RHS_var_indx)); 
    for regr = 2:num_regrs_base; fprintf(p.fid, '%8.2f', TS.SD_SD_x_adj(1:2, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end
    for regr = num_regrs_base + 1:num_regrs; fprintf(p.fid, '%8.2f', TS.SD_SD_x_adj(1:3, regr, LHS_per_indx, LHS_RHS_var_indx)); fprintf(p.fid, '    '); end

    ave_resid_cor_shrt = squeeze(TS.ave_resid_cor(:, LHS_per_indx, LHS_RHS_var_indx));
    fprintf(p.fid, '\nAve resid cor%6.2f     ', ave_resid_cor_shrt(1)); 
    for regr = 2:num_regrs_base;  fprintf(p.fid, '%8.2f            ', ave_resid_cor_shrt(regr)); end
    for regr = num_regrs_base + 1:num_regrs;  fprintf(p.fid, '%8.2f                    ', ave_resid_cor_shrt(regr)); end

  end   % LHS_per
end

% OK
%
function TS_tmp = run_set_of_TS_dR_on_PR_and_one_dP(LHS_per, LHS_var, RHS_var, style_passed, c, d, p)

  % t = month
  % c = area
  %
  % K  = number of coefs, including intercept if regression has one
  % N  = total number of observations across metro areas and months
  % Nt = number of observations (metro areas with data) for month t
  % Nc = number of observations (months) for area i
  %
  %                                  
  % c.not_demeaned                   TS       SE2c = {sum t e2ct} / (Nc - K)                        SE2 = sum c {Nc SE2c / N}
  %								  
  % c.demean_CS_w_int or _wo_int     TS       SE2c = {sum t [Nt/(Nt - 1)]e2ct} / (Nc - K)           SE2 = sum c {Nc SE2c / N}
  %								  
  % 
  % frst_mth is frst form mth for area   num_regrs x 1
  % last_mth is last form mth for area   scalar
  % nt is number of areas in month   num_regrs x c.num_mths_ttl

  if ~x_in_y(style_passed, [c.not_demeaned, c.demean_CS_w_int, c.demean_CS_wo_int, c.demean_both_full])
    disp('run_set_of_TS_dR_on_PR_and_one_dP cannot handle this style_passed')
    c.demean_label(style_passed)
    pause
  end

  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.FM, style_passed, c, d, p);

  [num_coef_stats, num_area_coef_stats, num_ave_ts, num_area_ts] = get_num_TS_coef_stats(p);   % num_coef_stats is for the CS ave of multiple TS regressions, num_area_coef_stats is for one TS regression.


  % Global (cross-area) results
  %
  coef_stats             = NaN(num_coef_stats, max_num_coefs,  num_regrs);
  ave_SD_x_adj           =             NaN(max_num_slopes, num_regrs);
  SD_SD_x_adj            =             NaN(max_num_slopes, num_regrs);

  glbl_r2_ttl_raw        =                             NaN(num_regrs, 1);    % 1 - variance of bias adjusted residuals relative to variance of raw y = 1 - var_e/var_raw_y)
  glbl_r2_dmn_raw        =                             NaN(num_regrs, 1);    % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
  glbl_r2_slps_raw_y     =                             NaN(num_regrs, 1);    % fraction of raw variance explained by x = R2 total - R2 demean = 
  glbl_r2_slps_dmn_y     =                             NaN(num_regrs, 1);    % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y

  glbl_var_e             =                           zeros(num_regrs, 1);
  glbl_var_raw_y         =                           zeros(num_regrs, 1);
  glbl_var_BA_adj_y      =                           zeros(num_regrs, 1);
  
  glbl_cov_pnl_wght_coef = NaN(max_num_coefs, num_regrs);
  glbl_cor_pnl_wght_coef = NaN(max_num_coefs, num_regrs);
  glbl_cor_var_x_coef    = NaN(max_num_coefs, num_regrs);
  
  ave_r2_raw_y           =                             NaN(num_regrs, 1);
  ave_r2_dmn_y           =                             NaN(num_regrs, 1);
  ave_var_e              =                             NaN(num_regrs, 1);
  ave_resid_cor          =                             NaN(num_regrs, 1);
  se_resid_cor_good      =                             NaN(num_regrs, 1);

  % Area results
  %
  coef_stats_area        = NaN(num_area_coef_stats, max_num_coefs, num_regrs, c.num_areas_both);
  autos_area             =                          NaN(num_autos, num_regrs, c.num_areas_both);
  r2_raw_y_area          =                                     NaN(num_regrs, c.num_areas_both);
  r2_dmn_y_area          =                                     NaN(num_regrs, c.num_areas_both);
  var_e_area             =                                     NaN(num_regrs, c.num_areas_both);
  SD_x_adj_area          =                     NaN(max_num_slopes, num_regrs, c.num_areas_both);
  cor_PR_dP_area         =                                  NaN(num_RHS_pers, c.num_areas_both);
  

  % Working arrays
  %
  y_raw_ttl    = NaN(c.max_TS_CS_obs, 1);
  y_BA_adj_ttl = NaN(c.max_TS_CS_obs, 1);

  y_raw     = NaN(c.last_mth, 1);
  y_adj     = NaN(c.last_mth, 1);
  y_BA_adj  = NaN(c.last_mth, 1);

  x_raw     = NaN(c.last_mth, max_num_coefs);
  x_adj     = NaN(c.last_mth, max_num_coefs);
  resids_FE = NaN(c.last_mth, 1);

  sum_x2    = NaN(c.num_areas_both, max_num_coefs);
  var_x     = NaN(c.num_areas_both, max_num_coefs);
  coef_area = NaN(c.num_areas_both, max_num_coefs);


  if x_in_y(style_passed, [c.not_demeaned, c.demean_CS_w_int])
    x_adj(:, 1) = 1;
    x_raw(:, 1) = 1;
  end


  % determine frst and last form mth for each area and each regr
  %
  % frst_form_mth = NaN(num_regrs, 1);  
  % last_form_mth = NaN;                       
  % frst_form_mth_area = NaN(c.num_areas_both, num_regrs);    % 1 + num_RHS_pers...1 is for simple P/R
  % last_form_mth_area = NaN(c.num_areas_both, 1);            % 1 + num_RHS_pers
  %
  % nt = zeros(c.num_mths_ttl,    num_RHS_pers + 1)            % number of good areas each month
  % nc = zeros(c.num_areas_both, num_RHS_pers + 1)            % number of good months for each area
  %
  [frst_form_mth, last_form_mth, frst_form_mth_area, last_form_mth_area, nt, nc] = determine_and_wrt_frst_and_last_form_mth(LHS_per, LHS_var, RHS_var, c.dont_wrt_frst_last_mth, c, d, p);
  

  % Select LHS_chg data based on LHS_var and LHS_per...0 means for all areas
  %
  LHS_raw = load_mthly_qtrly_semi_or_annual_chg (LHS_var, LHS_per, 0, c, d, p);  % The two omitted results are for subtracting the overall mean for the month (in the TS regression)


  % Load vector of PR or RP data
  %
  PR_raw = load_full_PR_vector_regr(0, c, d, p);

  if style_passed == c.not_demeaned
    LHS_adj   = LHS_raw;
    PR_adj    = PR_raw;
    resids_FE = LHS_raw;
  end

  % Run TS regression with only PR if regr = 1 or PR and dP if regr > 1
  %
  for regr = 1:num_regrs

    % Setup RHS variable
    %
    if regr == 1
      last_slope = frst_slope;
      num_slopes = 1;
      regr_base  = regr;

    else
      last_slope = frst_slope + 1;
      num_slopes = 2;
      regr_base  = regr;
    end


    if regr ~= 1
      RHS_per      = p.RHS_per(regr_base - 1);
      RHS_raw_lag  = load_mthly_qtrly_semi_or_annual_chg(RHS_var, RHS_per, 0, c, d, p);  
      RHS_raw_lead = load_mthly_qtrly_semi_or_annual_chg(RHS_var, LHS_per, 0, c, d, p);  
      
      if style_passed == c.not_demeaned
        RHS_adj_lag = RHS_raw_lag;
        
      elseif style_passed == c.demean_both_full
        RHS_adj_lag =  demean_x_both_full(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, RHS_raw_lag, c, d, p);

      else
        RHS_adj_lag = demean_x_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, RHS_raw_lag, c, d, p);

      end
    end

    if style_passed == c.demean_both_full
      LHS_adj = demean_y_both_full(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, LHS_raw, c, d, p);
      PR_adj  = demean_x_both_full(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, PR_raw, c, d, p);
    
      % adjust resids_FE for CS demeaning bias
      %
      if regr == 1
        resids_FE = LHS_adj;
        for area = 1:c.num_areas_both
          for mth = frst_form_mth(1):last_form_mth
            resids_FE(mth, area) = resids_FE(mth, area) * sqrt(nt(mth, regr_base) / (nt(mth, regr_base) - 1)) * sqrt(nc(area, regr_base) / (nc(area, regr_base) - 1));
          end
        end
      end


    elseif style_passed ~= c.not_demeaned
      LHS_adj = demean_y_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, LHS_per, LHS_raw, c, d, p);
      PR_adj  = demean_x_each_month(frst_form_mth(regr_base), last_form_mth, frst_form_mth_area(:, regr_base), last_form_mth_area, PR_raw, c, d, p);
      
      % adjust resids_FE for CS demeaning bias
      %
      if regr == 1
        resids_FE = LHS_adj;
        for mth = frst_form_mth(1):last_form_mth
          resids_FE(mth, :) = resids_FE(mth, :) * sqrt(nt(mth, regr_base) / (nt(mth, regr_base) - 1));
        end
      end

    end

    resids_area = NaN(c.num_mths_ttl, c.num_areas_both);
    
    
    % Compute TS regressions for each area with only PR if regr = 1 or PR and dP if regr > 1
    %
    SSE         = 0;
    SSE_FE      = 0;
    num_obs_ttl = 0;

    for area = 1:c.num_areas_both
      frst_mth = frst_form_mth_area(area, regr_base);
      last_mth = last_form_mth_area(area);
      
      y_raw(frst_mth:last_mth) = LHS_raw(frst_mth + LHS_per:last_mth + LHS_per, area);
      y_adj(frst_mth:last_mth) = LHS_adj(frst_mth + LHS_per:last_mth + LHS_per, area);
    
      if style_passed == c.not_demeaned
        y_BA_adj(frst_mth:last_mth) = y_adj(frst_mth:last_mth);
        
      elseif style_passed == c.demean_both_full
        y_BA_adj(frst_mth:last_mth) = y_adj(frst_mth:last_mth) .* sqrt(nt(frst_mth:last_mth, regr_base) ./ (nt(frst_mth:last_mth, regr_base) - 1)) .* sqrt(nc(area, regr_base) / (nc(area, regr_base) - 1));

      else
        y_BA_adj(frst_mth:last_mth) = y_adj(frst_mth:last_mth) .* sqrt(nt(frst_mth:last_mth, regr_base) ./ (nt(frst_mth:last_mth, regr_base) - 1));
      end

      x_raw(frst_mth:last_mth, frst_slope) =  PR_raw(frst_mth:last_mth, area);
      x_adj(frst_mth:last_mth, frst_slope) =  PR_adj(frst_mth:last_mth, area);
      SD_x_adj_area(1, regr, area)         = std(x_adj(frst_mth:last_mth, frst_slope));

      
      if regr ~= 1
        x_raw(frst_mth:last_mth, frst_slope + 1) = RHS_raw_lag(frst_mth:last_mth, area);
        x_adj(frst_mth:last_mth, frst_slope + 1) = RHS_adj_lag(frst_mth:last_mth, area);

        SD_x_adj_area(2, regr, area) = std(x_adj(frst_mth:last_mth, frst_slope + 1));
      end
      
      if 1 < regr && regr <= 1 + num_RHS_pers
        cor_PR_dP_area(regr - 1, area) = corr(x_adj(frst_mth:last_mth, frst_slope), x_adj(frst_mth:last_mth, frst_slope + 1));   % num_RHS_pers x num_areas_both
      end


      % Run TS regression with only PR if regr = 1, PR and dP otherwise
      %
      [coef_stats_shrt, r2_raw_y_area_shrt, r2_dmn_y_area_shrt, var_e_shrt, resids_adj, autos_shrt] = TS_reg_w_autos(x_adj(frst_mth:last_mth, 1:last_slope), y_raw(frst_mth:last_mth), y_adj(frst_mth:last_mth), y_BA_adj(frst_mth:last_mth), LHS_per, style_passed, nt(frst_mth:last_mth, regr_base), c, d, p);


      % Store sum of x2 and var of x to look at relation between coefs and N*Var_x or Var_x  Note, Var_x is max likelihood value
      %
      if style_passed == c.demean_CS_w_int
        sum_x2(area, frst_slope:last_slope) = sum(x_adj(frst_mth:last_mth, frst_slope:last_slope).^2); 
         var_x(area, frst_slope:last_slope) = sum_x2(area, frst_slope:last_slope) ./ (last_mth - frst_mth - 1);
      end
        
  
      % compute implied intercept and store in coef_stats_area if regression does not include it
      %
      if style_passed == c.demean_CS_wo_int || style_passed == c.demean_both_full
        coef_stats_area(1, 1,                regr, area) = mean(y_raw(frst_mth:last_mth, 1)) - mean(x_raw(frst_mth:last_mth, 1:last_slope)) * coef_stats_shrt(1, 1:last_slope)';
        coef_stats_area(:, 2:last_slope + 1, regr, area) = coef_stats_shrt;
      else
        coef_stats_area(:, 1:last_slope, regr, area) = coef_stats_shrt;
      end
	
	  coef_area(area, 1:last_slope) = coef_stats_area(1, 1:last_slope, regr, area);
	   
      resids_area(frst_mth:last_mth, area) = resids_adj;
      
      autos_area(:, regr, area) = autos_shrt;
      r2_raw_y_area(regr, area) = r2_raw_y_area_shrt;
      r2_dmn_y_area(regr, area) = r2_dmn_y_area_shrt;
         var_e_area(regr, area) = var_e_shrt;
	
	
      % Cumulate SSE and store y_raw for all observations
      %
      SSE = SSE + (nc(area, regr_base) / (nc(area, regr_base) - last_slope)) * sum(resids_adj.^2, 1);
                    
         y_raw_ttl(num_obs_ttl + 1:num_obs_ttl + nc(area, regr_base)) =    y_raw(frst_mth:last_mth);
      y_BA_adj_ttl(num_obs_ttl + 1:num_obs_ttl + nc(area, regr_base)) = y_BA_adj(frst_mth:last_mth);
           
      num_obs_ttl = num_obs_ttl + nc(area, regr_base);

    end
    
    if style_passed == c.demean_CS_w_int
      pnl_wght(:, frst_slope:last_slope) = sum_x2(:, frst_slope:last_slope) ./ sum(sum_x2(:, frst_slope:last_slope));

      for slp = frst_slope:last_slope
        Tcov_pnl_wght_coef_tmp            = c.num_areas_both * cov(pnl_wght(:, slp), coef_area(:, slp));
        glbl_cov_pnl_wght_coef(slp, regr) = Tcov_pnl_wght_coef_tmp(1, 2);
        glbl_cor_pnl_wght_coef(slp, regr) = corr(pnl_wght(:, slp), coef_area(:, slp));
           glbl_cor_var_x_coef(slp, regr) =    corr(var_x(:, slp), coef_area(:, slp));
      end
    end    

        glbl_var_e(regr) = SSE / num_obs_ttl;
    glbl_var_raw_y(regr) = var(y_raw_ttl(1:num_obs_ttl));

    if style_passed == c.not_demeaned
      glbl_var_BA_adj_y(regr) = glbl_var_raw_y(regr);  
    else
      glbl_var_BA_adj_y(regr) = sum(y_BA_adj_ttl(1:num_obs_ttl).^2) ./ num_obs_ttl;  % demeaned monthly, so can sum squared values, adjusted for overfitting bias above divide by num_obs_ttl
    end

       glbl_r2_ttl_raw(regr) = 1 - glbl_var_e(regr)        / glbl_var_raw_y(regr);
       glbl_r2_dmn_raw(regr) = 1 - glbl_var_BA_adj_y(regr) / glbl_var_raw_y(regr);     % fraction of raw variance explained by demeaning = 1 - variance of bias adjusted demeaned y / variance of y = 1 - var_BA_adj_y / var_raw_y
    glbl_r2_slps_dmn_y(regr) = 1 - glbl_var_e(regr)        / glbl_var_BA_adj_y(regr);  % fraction of demeaned variance explained by x = 1 - variance of bias adjusted residuals relative to bias adjusted demeaned  y
    glbl_r2_slps_raw_y(regr) = glbl_r2_ttl_raw(regr) - glbl_r2_dmn_raw(regr);          % fraction of raw variance explained by x = R2 total - R2 demean

    ave_r2_raw_y(regr) = mean(r2_raw_y_area(regr, :));
    ave_r2_dmn_y(regr) = mean(r2_dmn_y_area(regr, :));
       ave_var_e(regr) = mean(var_e_area(regr, :));
    
    ave_SD_x_adj(1:num_slopes, regr) = mean(SD_x_adj_area(1:num_slopes, regr, :), 3, 'omitnan');
     SD_SD_x_adj(1:num_slopes, regr) =  std(SD_x_adj_area(1:num_slopes, regr, :), 0, 3, 'omitnan');
     
    [coef_stats(:, 1:num_slopes + 1, regr), ave_resid_cor(regr), se_resid_cor_good(regr)] = summarize_TS_area_coefs(squeeze(coef_stats_area(1, 1:num_slopes + 1, regr, :))', resids_area, c, d, p);
  end

  TS_tmp.coef_stats         = coef_stats;  
  TS_tmp.se_resid_cor_good  = se_resid_cor_good;  
  TS_tmp.ave_resid_cor      = ave_resid_cor;  
  TS_tmp.glbl_var_e         = glbl_var_e;
  TS_tmp.glbl_var_raw_y     = glbl_var_raw_y;
  TS_tmp.glbl_var_BA_adj_y  = glbl_var_BA_adj_y;
  TS_tmp.glbl_r2_ttl_raw    = glbl_r2_ttl_raw;
  TS_tmp.glbl_r2_dmn_raw    = glbl_r2_dmn_raw;
  TS_tmp.glbl_r2_slps_raw_y = glbl_r2_slps_raw_y;
  TS_tmp.glbl_r2_slps_dmn_y = glbl_r2_slps_dmn_y;
  TS_tmp.ave_r2_raw_y       = ave_r2_raw_y;
  TS_tmp.ave_r2_dmn_y       = ave_r2_dmn_y;
  TS_tmp.ave_var_e	        = ave_var_e;
  TS_tmp.ave_SD_x_adj       = ave_SD_x_adj;
  TS_tmp.SD_SD_x_adj        = SD_SD_x_adj;
  TS_tmp.cor_PR_dP          = mean(cor_PR_dP_area', 'omitnan')';  % cor_PR_dp_area is num_RHS_pers x num_areas_both,             ave_cor_PR_dP is num_RHS_pers x 1


  if style_passed == c.demean_CS_w_int
    TS_tmp.glbl_cov_pnl_wght_coef = glbl_cov_pnl_wght_coef;
    TS_tmp.glbl_cor_pnl_wght_coef = glbl_cor_pnl_wght_coef;
    TS_tmp.glbl_cor_var_x_coef    = glbl_cor_var_x_coef;
  else
    TS_tmp.glbl_cov_pnl_wght_coef = NaN;
    TS_tmp.glbl_cor_pnl_wght_coef = NaN;
    TS_tmp.glbl_cor_var_x_coef    = NaN;
  end
 
  TS_tmp.coef_stats_area   = coef_stats_area;
  TS_tmp.autos_area        = autos_area;
  TS_tmp.r2_raw_y_area     = r2_raw_y_area;
  TS_tmp.r2_dmn_y_area     = r2_dmn_y_area;
  TS_tmp.var_e_area        = var_e_area;
  TS_tmp.SD_x_adj_area     = SD_x_adj_area;
  TS_tmp.cor_PR_dP_area    = cor_PR_dP_area;           % cor_PR_dp_area is num_RHS_pers x num_areas_both
  
end


function [coef_stats, r2_raw_y_area, r2_dmn_y_area, var_e, resids_adj, autos] = TS_reg_w_autos(x_adj, y_raw, y_adj, y_BA_dmn, y_per, style_passed, nt, c, d, p)

  % y_raw and y_adj are num_obs x 1
  % y_raw is y before subtracting CS average...for computing R2
  % y_adj is y after subtracting CS average...for computing regr
  %
  % x is num_obs x num_coefs
  %
  % coef_stats is 7 x num_coefs:  coef, iid SE, hetero SE, area clstrd SE, t iid, t hetero, t area clustrd
  % r2 is scalar
  % autos is num_autos x 1 = p.num_autos_resids x 1
  %
  num_obs   = size(x_adj, 1);
  num_coefs = size(x_adj, 2);
  num_autos = max([p.num_autos_resids, p.num_autos_in_TS_coef_se, y_per - 1]);

  
  % Compute coefficients
  %
  coef = (x_adj' * x_adj) \ x_adj' * y_adj;
  %Kx1    KxT      TxK      KxT      Tx1
  
  coef = coef';
  %1xK
  

  % The CS average for each month may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y_raw = var(y_raw);
  var_y_dmn = var(y_adj);   % No DOF adjustment needed because mean of 

  if style_passed == c.demean_both_full   % mean already set to zero and BA adjusts for degrees of freedom used
    var_y_BA_dmn = sum(y_BA_dmn.^2) / num_obs;
  else                                    % either not demeaning or CS demeaned  var divides by n-1
    var_y_BA_dmn = var(y_BA_dmn);
  end  

  % Calculate residuals and adjust for overfitting if demeaned
  %
  resids = y_adj - x_adj * coef';
  %         Tx1     TxK     Kx1
           
%  if x_in_y(style_passed, [c.demean_CS_w_int, c.demean_CS_wo_int])
  if style_passed ~= c.not_demeaned
    resids_adj = sqrt(nt ./(nt - 1)) .* resids;                            % adjust for bias caused by CS demeaning
  else
    resids_adj = resids;
  end

  % Calculate residual variances and r2
  %
  var_e = sum(resids_adj.^2, 1) / (num_obs - num_coefs);
  r2_raw_y_area = 1.0 - (var_e ./ var_y_raw);          % 1 x 1
  r2_dmn_y_area = 1.0 - (var_e ./ var_y_dmn);          % 1 x 1


  autos_shrt = autocorr(resids_adj, 'NumLags', num_autos);
  autos      = autos_shrt(2:num_autos + 1, 1);    % num_autos x 1


  % Compute SEs of coefs assuming iid resids
  %
  xx_adj     = x_adj' * x_adj;                                     % num_coefs x num_coefs
  inv_xx_adj = inv(xx_adj);

  iid_cov_TS = var_e * inv_xx_adj;
  %   KxK       1x1      K x K


  % sum autos up to y_per - 1 
  %
  if y_per == 1
    auto_adj = 1;
  else
    num_lags = y_per - 1;
    auto_adj = cmpt_auto_adj(num_lags, num_obs, autos(1:num_lags));
  end

  
  % Heteroscedastic consistent SEs
  %
%  disp("Hetero")
%  z_tmp = diag(resids_adj.^2);
%  size(z_tmp)
%  z_tmp(1:10, 1:10)
  
  x_omega_x_white = x_adj' * diag(resids_adj.^2) * x_adj;            % diag makes a square matrix with resid_adj^2 on the diagonal

  white_cov_TS = inv_xx_adj * x_omega_x_white * inv_xx_adj;    % num_coefs x num_coefs  
  

  % Generalized HH
  %
  sig2 = resids_adj * resids_adj';
  
  for row = 1:num_obs
    for col = 1:row - y_per
      sig2(row, col) = 0;
      sig2(col, row) = 0;
    end
  end
  
%  disp("sig2")
%  y_per
%  size(sig2)
%  sig2(1:10, 1:10)
%  pause
  
  genrlzd_HH_x_sig2_x = x_adj' * sig2 * x_adj;
  %                      KxT     T x T   TxK
      
  genrlzd_HH_cov_TS = inv_xx_adj * genrlzd_HH_x_sig2_x * inv_xx_adj;
  
  if p.num_autos_in_TS_coef_se == 0
        SE_TS = NaN(4, num_coefs);
    t_stat_TS = NaN(4, num_coefs);
  else
        SE_TS = NaN(4 + size(p.num_autos_in_TS_coef_se, 2), num_coefs);
    t_stat_TS = NaN(4 + size(p.num_autos_in_TS_coef_se, 2), num_coefs);
  end
  
  SE_TS(1:4, :) = [diag(iid_cov_TS)'; auto_adj * diag(iid_cov_TS)'; diag(white_cov_TS)'; diag(genrlzd_HH_cov_TS)'];
  % 4xK
  
  if p.num_autos_in_TS_coef_se ~= 0
    for auto_indx = 1:size(p.num_autos_in_TS_coef_se, 2)
      num_lags = p.num_autos_in_TS_coef_se(auto_indx);
      auto_adj = cmpt_auto_adj(num_lags, num_obs, autos(1:num_lags));
      SE_TS(4 + auto_indx, :) = auto_adj * diag(iid_cov_TS)';
    end
  end
  
  SE_TS = sqrt(SE_TS);

  for indx = 1:size(SE_TS, 1)
    t_stat_TS(indx, :) = coef ./ SE_TS(indx, :);
  end
  
  coef_stats = [coef; SE_TS; t_stat_TS];         % coefs, se_coefs, t-stats  1 + 2 * (4 + size(p.num_autos_in_TS_coef_se, 2)) x num_coefs
  
end


function [TS_reg_coefs, TS_reg_r2, TS_reg_var_y, TS_reg_var_e] = TS_reg_smple(x, y, y_per, c, d, p)

  % y is num_obs x 1
  % x is num_obs x num_coefs
  %                                 1    2       3       4          5               6      7         8         9
  % coef_stats is 9 x num_coefs:  coef, iid SE, LHS-1, hetero SE, generalized HH, t iid, t LHS-1, t hetero, t generalized HH
  % r2 is scalar
  num_obs   = size(x, 1);
  num_coefs = size(x, 2);
  
  
  % Compute coefficients
  %
  coef = (x' * x) \ x' * y;
  %Kx1    KxT TxK   KxT Tx1
  
  coef = coef';
  %1xK
  

  % The CS average for each month may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y = var(y);


  % Calculate residuals and adjust for overfitting if demeaned
  %
  resids = y - x * coef';
  %       Tx1 TxK   Kx1
           

  % Calculate residual variances and r2
  %
  var_e = sum(resids.^2, 1) / (num_obs - num_coefs);
  r2    = 1.0 - (var_e ./ var_y);          % 1 x 1

  % Compute SEs of coefs assuming iid resids
  %
  xx     = x' * x;                                     % num_coefs x num_coefs
  inv_xx = inv(xx);

  iid_cov_TS = var_e * inv_xx;
  %   KxK       1x1    K x K


  % sum autos up to y_per - 1 
  %
  if y_per == 1
    auto_adj = 1;
  else
    num_lags   = y_per - 1;
    autos_shrt = autocorr(resids, 'NumLags', num_lags);
    autos      = autos_shrt(2:num_lags + 1, 1);    % num_autos x 1

    auto_adj   = cmpt_auto_adj(num_lags, num_obs, autos(1:num_lags));
  end

  
  % Heteroscedastic consistent SEs
  %
%  disp("Hetero")
%  z_tmp = diag(resids_adj.^2);
%  size(z_tmp)
%  z_tmp(1:10, 1:10)
  
  x_omega_x_white = x' * diag(resids.^2) * x;          % diag makes a square matrix with resid_adj^2 on the diagonal

  white_cov_TS = inv_xx * x_omega_x_white * inv_xx;    % num_coefs x num_coefs  
  

  % Generalized HH
  %
  sig2 = resids * resids';
  
  for row = 1:num_obs
    for col = 1:row - y_per
      sig2(row, col) = 0;
      sig2(col, row) = 0;
    end
  end
  
  genrlzd_HH_x_sig2_x = x' * sig2 * x;
  %                    KxT  T x T  TxK
      
  genrlzd_HH_cov_TS = inv_xx * genrlzd_HH_x_sig2_x * inv_xx;
  
      SE_TS = NaN(4, num_coefs);
  t_stat_TS = NaN(4, num_coefs);
  
  SE_TS(1:4, :) = [diag(iid_cov_TS)'; auto_adj * diag(iid_cov_TS)'; diag(white_cov_TS)'; diag(genrlzd_HH_cov_TS)'];
  % 4xK
  
  SE_TS = sqrt(SE_TS);

  for indx = 1:size(SE_TS, 1)
    t_stat_TS(indx, :) = coef ./ SE_TS(indx, :);
  end
  
  TS_reg_coefs = NaN(9, 4);
  TS_reg_coefs(:, 1:num_coefs) = [coef; SE_TS; t_stat_TS];         % coefs, se_coefs, t-stats  1 + 4 + 4
  TS_reg_r2 = r2;
  TS_reg_var_e = var_e;
  TS_reg_var_y = var_y;
  
end

% Needed
%
function [TS_reg_coefs, TS_reg_r2, TS_reg_var_y, TS_reg_var_e] = TS_reg_FEs(x, y, y_per, cnt, LHS_obs, c, d, p)

  % y is num_obs x 1
  % x is num_obs x num_coefs
  %                                 1    2       3       4          5               6      7         8         9
  % coef_stats is 9 x num_coefs:  coef, iid SE, LHS-1, hetero SE, generalized HH, t iid, t LHS-1, t hetero, t generalized HH
  % r2 is scalar
  num_obs   = size(x, 1);
  num_coefs = size(x, 2);
  ttl_obs   = size(LHS_obs, 1);

  
  
  % Compute coefficients
  %
  coef = (x' * x) \ x' * y;
  %Kx1    KxT TxK   KxT Tx1
  
  coef = coef';
  %1xK
  

  % The CS average for each month may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y = var(y);


  % Calculate residuals and adjust for overfitting if demeaned
  %
  resids = y - x * coef';
  %       Tx1 TxK   Kx1
           

  % Calculate residual variances and r2   cnt is the number of observations in a month or area
  %
  sse = 0;
  for grp = 1:num_obs
    sse = sse + cnt(grp) * resids(grp).^2;
  end
    
  var_e = (sse / ttl_obs) * (num_obs / (num_obs - num_coefs));   % Adjust for overfitting bias...num_areas / (num_areas - 1) or num_mths / (num_mths - 1)
  
  r2    = 1.0 - (var_e ./ var_y);          % 1 x 1

  % Compute SEs of coefs assuming iid resids
  %
  xx     = x' * x;                                     % num_coefs x num_coefs
  inv_xx = inv(xx);

  iid_cov_TS = var_e * inv_xx;
  %   KxK       1x1    K x K


  % sum autos up to y_per - 1 
  %
  if y_per == 1
    auto_adj = 1;
  else
    num_lags   = y_per - 1;
    autos_shrt = autocorr(resids, 'NumLags', num_lags);
    autos      = autos_shrt(2:num_lags + 1, 1);    % num_autos x 1

    auto_adj   = cmpt_auto_adj(num_lags, num_obs, autos(1:num_lags));
  end

  
  % Heteroscedastic consistent SEs
  %
%  disp("Hetero")
%  z_tmp = diag(resids_adj.^2);
%  size(z_tmp)
%  z_tmp(1:10, 1:10)
  
  x_omega_x_white = x' * diag(resids.^2) * x;          % diag makes a square matrix with resid_adj^2 on the diagonal

  white_cov_TS = inv_xx * x_omega_x_white * inv_xx;    % num_coefs x num_coefs  
  

  % Generalized HH
  %
  sig2 = resids * resids';
  
  for row = 1:num_obs
    for col = 1:row - y_per
      sig2(row, col) = 0;
      sig2(col, row) = 0;
    end
  end
  
  genrlzd_HH_x_sig2_x = x' * sig2 * x;
  %                    KxT  T x T  TxK
      
  genrlzd_HH_cov_TS = inv_xx * genrlzd_HH_x_sig2_x * inv_xx;
  
      SE_TS = NaN(4, num_coefs);
  t_stat_TS = NaN(4, num_coefs);
  
  SE_TS(1:4, :) = [diag(iid_cov_TS)'; auto_adj * diag(iid_cov_TS)'; diag(white_cov_TS)'; diag(genrlzd_HH_cov_TS)'];
  % 4xK
  
  SE_TS = sqrt(SE_TS);

  for indx = 1:size(SE_TS, 1)
    t_stat_TS(indx, :) = coef ./ SE_TS(indx, :);
  end
  
  TS_reg_coefs = NaN(9, 3);
  TS_reg_coefs(:, 1:num_coefs) = [coef; SE_TS; t_stat_TS];         % coefs, se_coefs, t-stats  1 + 4 + 4
  TS_reg_r2 = r2;
  TS_reg_var_e = var_e;
  TS_reg_var_y = var_y;
  
end


function [CS_reg_coefs, CS_reg_r2, CS_reg_var_y, CS_reg_var_e] = CS_reg_smple(x, y, c, d, p)
           % 3 X K         1 X 1
           
  num_obs   = size(x, 1);
  num_coefs = size(x, 2);
  
  
  % Compute coefficients
  %
  coef = (x' * x) \ x' * y;
  %Kx1    KxT TxK   KxT Tx1
  
  coef = coef';
  %1xK
  

  % The CS average for each month may be subtracted in y_adj. Calculate variance of unadjusted y
  %
  var_y = var(y);


  % Calculate residuals and adjust for overfitting if demeaned
  %
  resids = y - x * coef';
  %       Tx1 TxK   Kx1
           

  % Calculate residual variances and r2
  %
  var_e = sum(resids.^2, 1) / (num_obs - num_coefs);
  r2    = 1.0 - (var_e ./ var_y);          % 1 x 1

  % Compute SEs of coefs assuming iid resids
  %
  xx     = x' * x;                                     % num_coefs x num_coefs
  inv_xx = inv(xx);

  iid_cov_TS = var_e * inv_xx;
  %   KxK       1x1    K x K


  SE_CS     = NaN(1, num_coefs);
  t_stat_CS = NaN(1, num_coefs);
  
  SE_CS = sqrt(diag(iid_cov_TS)');
  % 1xK
  
  t_stat_CS = coef ./ SE_CS;
  
  CS_reg_coefs = NaN(3, 3);
  CS_reg_coefs(:, 1:num_coefs) = [coef; SE_CS; t_stat_CS];         % coefs, se_coefs, t-stats  1 + 1 + 1
  CS_reg_r2    = r2;
  CS_reg_var_y = var_y;
  CS_reg_var_e = var_e;
  
end


%
% *****************************************************************************************
%                                                                                         *
%                                   End TS Regressions                                    *
%                                                                                         *
% *****************************************************************************************
%
% *****************************************************************************************
%                                                                                         *
%                                   Write Summary Table                                   *
%                                                                                         *
% *****************************************************************************************
% Needed
%
function control_wrt_summary_tbl_202212(FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS, c, d, p)

  fprintf(p.fid, '\n\n\nPnl One          Panel regression, not demeaned, with one intercept');
  fprintf(p.fid, '\nTS Area          Time series regressions, not demeaned, with area intercepts');
  fprintf(p.fid, '\nPnl Mth          Panel regression, demeaning each month by its CS average of raw X and Y, regression does not have an intercept');
  fprintf(p.fid, '\nTS Mth&Area      Time series regressions, demeaning each month by its CS average of raw X and Y, with area intercepts');
  fprintf(p.fid, '\nFM Mth           FM cross-section regressions, not demeaned, with monthly intercepts');
  fprintf(p.fid, '\nFM DM wo         FM cross-section regressions, CS demeaned, without monthly intercepts');
  fprintf(p.fid, '\n');

  for LHS_RHS_var_indx = 1:size(p.LHS_RHS_vars, 2)
    wrt_summary_tbl_202212_header_dR_dPK(LHS_RHS_var_indx, c, p)

    for LHS_per_indx = [2, 4]  % 1:size(p.LHS_per, 2)
      wrt_summary_table_202212_dR_dPK(LHS_per_indx, LHS_RHS_var_indx, FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS, c, d, p);
    end
  end
  
end
      
% Needed
%
function wrt_summary_tbl_202212_header_dR_dPK(LHS_RHS_var_indx, c, p)

  LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
  RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);

  fprintf(p.fid, '\n\nRegressing %s on %s and %s_Lag   (wrt_summary_tbl_202212)', strtrim(load_dP_or_dR_label(LHS_var, c)), strtrim(load_PR_label_shrt(c, p)), load_dP_or_dR_label(RHS_var, c));
  fprintf(p.fid, '\n\nTS Both wo and FM DM CS wo are included because they give the correct split of R2 Ttl between Dmn and X.');
  fprintf(p.fid, '\n');      
  fprintf(p.fid, '\n          '); for RHS_per = p.summary_RHS_pers; fprintf(p.fid, '                                   RHS Period = %2d                          ', RHS_per); end
  fprintf(p.fid, '\n          '); for RHS_per = p.summary_RHS_pers; fprintf(p.fid, '              Slope            t-stat        Ttl R2   Ttl R2 Split      R2 X'); end
  fprintf(p.fid, '\n          '); for RHS_per = p.summary_RHS_pers; fprintf(p.fid, '           PR      dP        PR      dP      Raw Y     Dmn      X      Dmn Y'); end
%                                                                                                                            123456789012345678123456781234567890
end                                                                              

% Needed
%                                                                                                   
function wrt_summary_table_202212_dR_dPK(LHS_per_indx, LHS_RHS_var_indx, FM_not_DM, FM_DM_CS_wo, Panel_not_DM, Panel_DM_CS_wo_int, TS_not_DM, TS_DM_CS_w, FE_on_FE_TS, c, d, p)
  LHS_var = p.LHS_RHS_vars(1,   LHS_RHS_var_indx);
  RHS_var = p.LHS_RHS_vars(2:3, LHS_RHS_var_indx);

  fprintf(p.fid, '\n\nLHS = %s(t,t+%d)', load_dP_or_dR_label(LHS_var, c), p.LHS_per(LHS_per_indx));

  wrt_summary_line_dR_dPK(Panel_not_DM,       1, 9, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % Panel One        2 = Int           9 = Generalized HH
  wrt_summary_line_dR_dPK(TS_not_DM,          2, 5, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % TS Area          2 = W Int         7 = Resid CS cor (or iid if smaller) for TS slopes
  wrt_summary_line_dR_dPK(Panel_DM_CS_wo_int, 1, 9, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % Panel Mth        1 = No int        9 = Generalized HH
  wrt_summary_line_dR_dPK(TS_DM_CS_w,         2, 5, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % TS Mth w         2 = W Int         7 = Resid CS cor (or iid if smaller) for TS slopes
  wrt_summary_line_dR_dPK(FM_not_DM,          2, 5, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % FM Mth           2 = W Int         5 = LHS - 1 (or iid if smaller) for FM slopes
  wrt_summary_line_dR_dPK(FM_DM_CS_wo,        2, 5, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % FM CS DM wo      2 = W Int         5 = LHS - 1 (or iid if smaller) for FM slopes

  wrt_summary_line_dR_dPK(FE_on_FE_TS, 2, 9, LHS_per_indx, LHS_RHS_var_indx, c, d, p);    % FE on FE TS      2 = W Int

end

% Needed
%
function wrt_summary_line_dR_dPK(regr, frst_slp, t_stat, LHS_per_indx, LHS_RHS_var_indx, c, d, p)
  
  num_pers = size(p.summary_RHS_per_loc, 2);

  fprintf(p.fid, '\n%s', regr.label_4);
  for RHS_per_indx = 1:size(p.summary_RHS_pers, 2)
    RHS_loc = p.summary_RHS_per_loc(RHS_per_indx) + 1;

    fprintf(p.fid, '%10.2f%8.2f%10.2f%8.2f%10.2f%8.2f%8.2f%10.2f', regr.coef_stats(1, frst_slp:frst_slp + 1, RHS_loc, LHS_per_indx, LHS_RHS_var_indx), regr.coef_stats(t_stat, frst_slp:frst_slp + 1, RHS_loc, LHS_per_indx, LHS_RHS_var_indx), ...
    regr.glbl_r2_ttl_raw(RHS_loc, LHS_per_indx, LHS_RHS_var_indx), regr.glbl_r2_dmn_raw(RHS_loc, LHS_per_indx, LHS_RHS_var_indx), regr.glbl_r2_slps_raw_y(RHS_loc, LHS_per_indx, LHS_RHS_var_indx), regr.glbl_r2_slps_dmn_y(RHS_loc, LHS_per_indx, LHS_RHS_var_indx));

    fprintf(p.fid, '    '); 
  end
end
  

% Needed
%
function [FM_not_DM_out, FM_DM_CS_wo_out, Pnl_not_DM_out, Pnl_DM_CS_wo_out, TS_not_DM_out, TS_DM_CS_w_out, FE_on_FE_TS_out] = add_regr_data_labels(...
          FM_not_DM_in,  FM_DM_CS_wo_in,  Pnl_not_DM_in,  Pnl_DM_CS_wo_in,  TS_not_DM_in,  TS_DM_CS_w_in,  FE_on_FE_TS_in)

  FM_not_DM_out            = FM_not_DM_in;
  FM_DM_CS_wo_out          = FM_DM_CS_wo_in;
  Pnl_not_DM_out           = Pnl_not_DM_in;
  Pnl_DM_CS_wo_out         = Pnl_DM_CS_wo_in;
  TS_not_DM_out            = TS_not_DM_in;
  TS_DM_CS_w_out           = TS_DM_CS_w_in;
  FE_on_FE_TS_out          = FE_on_FE_TS_in;

  FM_not_DM_out.label      = 'FM, month intcept, nothing demeaned   ';
  FM_DM_CS_wo_out.label    = 'FM, w/o intcept, months demeaned      ';
  Pnl_not_DM_out.label     = 'Panel, intcept, nothing demeaned      ';
  Pnl_DM_CS_wo_out.label   = 'Panel, w/o intcept, months demeaned   ';
  TS_not_DM_out.label      = 'TS, area intcept, nothing demeaned    ';
  TS_DM_CS_w_out.label     = 'TS, area intcept, months demeaned     ';
  FE_on_FE_TS_out.label    = 'Mthly CS means regr on mthly CS means ';

  FM_not_DM_out.label_2    = 'FM      month     no  ';
  FM_DM_CS_wo_out.label_2  = 'FM       no     months';
  Pnl_not_DM_out.label_2   = 'Panel    yes      no  ';
  Pnl_DM_CS_wo_out.label_2 = 'Panel    no     months';
  TS_not_DM_out.label_2    = 'TS      area      no  ';
  TS_DM_CS_w_out.label_2   = 'TS      area    months';
  FE_on_FE_TS_out.label_2  = 'FE on FE TS regr      ';
												 
  FM_not_DM_out.label_3    = 'FM Mth  ';
  FM_DM_CS_wo_out.label_3  = 'FM CS wo';
  Pnl_not_DM_out.label_3   = 'Pnl One ';
  Pnl_DM_CS_wo_out.label_3 = 'Pnl CSwo';
  TS_not_DM_out.label_3    = 'TS Area ';
  TS_DM_CS_w_out.label_3   = 'TS CSw  ';
  FE_on_FE_TS_out.label_3  = 'FE TS   ';

  FM_not_DM_out.label_4    = 'FM Mth        ';
  FM_DM_CS_wo_out.label_4  = 'FM CS wo      ';
  Pnl_not_DM_out.label_4   = 'Panel One     ';
  Pnl_DM_CS_wo_out.label_4 = 'Pnl DM CS wo  ';
  TS_not_DM_out.label_4    = 'TS Area       ';
  TS_DM_CS_w_out.label_4   = 'TS CS w       ';
  FE_on_FE_TS_out.label_4  = 'TS on CS Means';
end

% Needed
%
function dum = cmpt_dum_r2(y, c, d, p)

  num_mths  = size(y, 1);
  num_areas = size(y, 2);
  sum_obs   = sum(~isnan(y), 'all');
  nc        = sum(~isnan(y), 1);  % number of good mths for each area (c is for city)
  nt        = sum(~isnan(y), 2);  % number of good areas for each mth (t is for time)

  y_TS_dmnd   = demean_y_TS_smpl(y, c, d, p);         
  y_CS_dmnd   = demean_y_CS_smpl(y, c, d, p);         
  y_both_dmnd = demean_y_both_smpl(y, c, d, p);
  
  var_y = var(y, 0, 'all', 'omitnan');
  
  sse = 0;
  for area = 1:num_areas
    sse = sse + (nc(area) / (nc(area) - 1)) * sum(y_TS_dmnd(:, area).^2, 'omitnan');
  end

  dum.r2_TS_dmnd = 1 - (sse / sum_obs) / var_y;
     
  sse = 0;
  for mth = 1:num_mths
    sse = sse + (nt(mth) / (nt(mth) - 1)) * sum(y_CS_dmnd(mth, :).^2, 'omitnan');
  end

  dum.r2_CS_dmnd = 1 - (sse / sum_obs) / var_y;
     
  sse = 0;
  for area = 1:num_areas
    for mth = 1:num_mths
      if ~isnan(y_both_dmnd(mth, area))
        sse = sse + (nt(mth) / (nt(mth) - 1)) * (nc(area) / (nc(area) - 1)) * y_both_dmnd(mth, area).^2;
      end
    end
  end

  dum.r2_both_dmnd = 1 - (sse / sum_obs) / var_y;
  
end     


% Needed
%
function FE_on_FE_TS = run_FE_on_FE_TS_regs(c, d, p)

  % ^^ Tables 1 and 3 ^^
  %
  % This function produces the TS regressions of CS Means on CS Means in Tables 1 and 3
  %
  [num_RHS_pers, max_RHS_per, num_LHS_pers, num_LHS_RHS_combos, num_regrs, num_regrs_base, num_autos, max_num_slopes, frst_slope, max_num_coefs] = setup_wrt_and_run_regr_constants(c.FE_TS, c.not_demeaned, c, d, p);
  num_coef_stats = get_num_FE_TS_coef_stats(p);                  

  FE_on_FE_TS.coef_stats      = NaN(num_coef_stats, max_num_coefs,  num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FE_on_FE_TS.glbl_r2_ttl_raw =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  %^^
  FE_on_FE_TS.glbl_r2_dmn_raw =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FE_on_FE_TS.glbl_r2_slps_raw_y =                              NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FE_on_FE_TS.glbl_r2_slps_dmn_y =                              NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);

  FE_on_FE_TS.glbl_var_e      =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FE_on_FE_TS.glbl_var_y      =                                 NaN(num_regrs, num_LHS_pers, num_LHS_RHS_combos);
  FE_on_FE_TS.cor_PR_dP       = NaN(num_RHS_pers, num_LHS_pers);

  max_ttl_obs = c.last_mth * c.num_areas_both;

   PR_CS_ave = NaN(c.last_mth, 1);
  LHS_CS_ave = NaN(c.last_mth, 1);
  RHS_CS_ave = NaN(c.last_mth, 1);
    cnt_area = NaN(c.last_mth, 1);


  % Cycle over LHS_RHS_combos and LHS periods
  %
  PR_raw  = load_full_PR_vector_regr(0, c, d, p);

  for LHS_RHS_var_indx = 1:num_LHS_RHS_combos
    LHS_var = p.LHS_RHS_vars(1, LHS_RHS_var_indx);
    RHS_var = p.LHS_RHS_vars(3, LHS_RHS_var_indx);

    for LHS_per_indx = 1:num_LHS_pers
      LHS_per = p.LHS_per(LHS_per_indx);

      LHS_raw = load_mthly_qtrly_semi_or_annual_chg(LHS_var, LHS_per, 0, c, d, p);

      for regr = 1:num_regrs

        x = NaN(p.last_mth, max_num_coefs);
        y = NaN(p.last_mth, 1);
  

        % Setup RHS variable
        %
        if regr == 1
          last_slope = frst_slope;
          num_slopes = 1;
          regr_base  = regr;

        else
          last_slope = frst_slope + 1;
          num_slopes = 2;
          regr_base  = regr;
        end

        ttl_obs = 0;
        LHS_obs = NaN(max_ttl_obs, 1);
        
        last_mth = p.last_mth - LHS_per;
          
        if regr == 1
          frst_mth = max(1 + p.RHS_per(1), 1 + p.PR_lag);
          for mth = frst_mth:last_mth
            LHS_CS_ave(mth) = mean(LHS_raw(mth + LHS_per, :), 'omitnan');
             PR_CS_ave(mth) =  mean(PR_raw(mth, :), 'omitnan');                                
               cnt_mth(mth) = sum(~isnan(PR_raw(mth, :)));
            for area = 1:c.num_areas_both
              if ~isnan(PR_raw(mth, area))
                ttl_obs = ttl_obs + 1;
                LHS_obs(ttl_obs) = LHS_raw(mth + LHS_per, area);
              end
            end
          end

        else
          RHS_per_indx = regr_base - 1;
          RHS_per      = p.RHS_per(RHS_per_indx);
          RHS_raw      = load_mthly_qtrly_semi_or_annual_chg(RHS_var, RHS_per, 0, c, d, p);  % This loads RHS_var for all areas

          frst_mth = max(1 + RHS_per, 1 + p.PR_lag);
            
          for mth = frst_mth:last_mth
            LHS_CS_ave(mth) = mean(LHS_raw(mth + LHS_per, :), 'omitnan');
            RHS_CS_ave(mth) = mean(RHS_raw(mth, :), 'omitnan');                                
             PR_CS_ave(mth) =  mean(PR_raw(mth, :), 'omitnan');                                
               cnt_mth(mth) = sum(~isnan(PR_raw(mth, :)));
            for area = 1:c.num_areas_both
              if ~isnan(PR_raw(mth, area))
                ttl_obs = ttl_obs + 1;
                LHS_obs(ttl_obs) = LHS_raw(mth + LHS_per, area);
              end
            end
          end
        end

        
        if regr == 1
          frst_mth = max(1 + p.RHS_per(1), 1 + p.PR_lag);
          obs      = last_mth - frst_mth + 1;

          [FE_on_FE_TS.coef_stats(:, :, 1, LHS_per_indx, LHS_RHS_var_indx), FE_on_FE_TS.glbl_r2_ttl_raw(1, LHS_per_indx, LHS_RHS_var_indx), FE_on_FE_TS.glbl_var_y(1, LHS_per_indx, LHS_RHS_var_indx), FE_on_FE_TS.glbl_var_e(1, LHS_per_indx, LHS_RHS_var_indx)] ...
           = TS_reg_FEs([ones(obs, 1), PR_CS_ave(frst_mth:last_mth)], LHS_CS_ave(frst_mth:last_mth), LHS_per, cnt_mth(frst_mth:last_mth), LHS_obs(1:ttl_obs), c, d, p);

        else
          frst_mth = max(1 + RHS_per, 1 + p.PR_lag);
          obs      = last_mth - frst_mth + 1;

          [FE_on_FE_TS.coef_stats(:, :, regr, LHS_per_indx, LHS_RHS_var_indx), FE_on_FE_TS.glbl_r2_ttl_raw(regr, LHS_per_indx, LHS_RHS_var_indx),  ...
                 FE_on_FE_TS.glbl_var_y(regr, LHS_per_indx, LHS_RHS_var_indx), FE_on_FE_TS.glbl_var_e(regr, LHS_per_indx, LHS_RHS_var_indx)]  = TS_reg_FEs([ones(obs, 1), PR_CS_ave(frst_mth:last_mth), RHS_CS_ave(frst_mth:last_mth)], LHS_CS_ave(frst_mth:last_mth), LHS_per, cnt_mth(frst_mth:last_mth), LHS_obs(1:ttl_obs), c, d, p);

          FE_on_FE_TS.cor_PR_dP(RHS_per_indx, LHS_per_indx) = corr(PR_CS_ave(frst_mth:last_mth), RHS_CS_ave(frst_mth:last_mth));
        end
      end
    end
  end
end



